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SUMMARY 


With the usual linearized supersonic theory, wings of given 
plan form and zero thickness are derived which have minimum 
drag for prescribed values of lift and pitching moment. The 
required wing twist and camber are represented as a summation 


of ten quasi-conical expressions, and the coefficients of this series 
The har- 


are determined by the usual calculus of variations 
monic method of Nikolskv, Kogan, Germain, and Graham is then 
used to determine the theoretically possible drag reduction to 


check the adequacy of the ten term series. 


Plan forms analyzed by the two methods are the delta and 


sweptback wings with both subsonic and supersonic edges. Re- 


sults indicate very large gains in (L/D),,,, for specific plan 


forms 


INTRODUCTION 


A‘ INTINUOUS TASK Of the applied aerodynamicist is 
to find a means of increasing the lift-drag ratio 


of high-performance aircraft and missiles, a task which 


is becoming increasingly difficult. The present paper 
will be confined to one facet of this problem and will 
consider the problem of optimizing the wing from the 
point of view of drag due to lift, which in most super- 


sonic aircraft forms an appreciable percentage of the 


total airplane drag. 

The present investigation will be based on the usual 
linearized theory and will be restricted to planar wing 
systems. Specifically we consider the problem of find- 


ing the generalized wing twist and camber which will 


yield minimum drag due to lift for a given wing plan 
form, flight Mach Number, lift, pitching moment, and 
possibly geometrical constraints on the wing surface, 


such as requiring straight control surface hinge lines. 


No attempt will be made here to find optimum plan- 
form shapes, although the harmonic method considered 
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here may be used to do so in an indirect manner for 
specific families of plan forms. 

The direct problem of finding the minimum drag due 
to lift and the determination of the necessary wing sur- 
face shape is not a new one (see, for example, references 
1 to 6). The procedure is straightforward, and it in- 
volves the consideration of a set of basic angle of attack 
distributions (e.g., the set composed of the constant 
angle of attack, linear twist, linear camber, etc.), the 
determination of their corresponding pressure distribu- 
tions, and then the finding of a suitable superposition 
of these basic distributions so as to obtain the desired 
optimum wing shape and drag by the usual minimi- 
zation procedure. The accessory requirements that the 
wing yields a given lift and pitching moment are satis- 
fied by using the Lagrangian multiplier method, while 
geometrical constraints on the wing surface shape are 
fulfilled by additional boundary conditions. The differ- 
ent papers referred to above differ primarily in the num- 
ber and generality of the basic distributions used. 

Although the above procedure reduces the optimi- 
zation procedure to ordinary calculus, it does introduce 
the question of the proper selection of a natural set of 
basic angle of attack distributions for a particular plan 
form, so that only a reasonable number of the set 1s 
required to represent the optimum wing shape. 

Supplementing the above procedure is a method of 
determining the actual minimum drag due to lift inde- 
pendent of any set of basic angle of attack distributions. 
This method, which we shall call the harmonic method, 
was first indicated by Nikolsky’ and given in detail by 
Kogan’ and Zhilin” using the control surface approach. 
Graham’ independently obtained the same harmonic 
problem by using the combined flow criterion for opti- 
mum drag obtained earlier by Jones." Ward'' and 
Germain,’ respectively, rederived and extended the 
above results to include linear constraints such as the 
pitching moment in addition to the lift. The harmonic 
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Fic. 1. Orientation of the coordinate system. 


method yields directly the minimum drag through the 
solution of a two-dimensional potential problem; but it 
is not, however, a straightforward calculation to find 
the required wing surface shape because it entails the 
solution of a nonclassical characteristic initial value 
problem. Graham** and Germain” have developed 
methods to compute the optimum wing shapes using 
this approach. Beane*! has applied Graham's method 
to a sonic edge diamond plan form. 

In the present paper no new approach is evolved. 
The somewhat inaccessible results of Nikolsky and 
Kogan and the brief note of Germain are reproduced 
and given in more detail; also, some recent results ob- 
tained by the Theoretical Aerodynamics Group, 
Convair (San Diego), are given to show typical gains in 
drag due to lift for several illustrative examples. 


DETERMINATION OF THE OPTIMUM WING SHAPE 


In the present paper we are concerned with those 
flows for which the Prandtl-Glauert equation and 
planar-boundary conditions are applicable. 

The basic flow equation is given by the wave equation 


—(M,? Ll) err + = 0 (1) 


where ¢ is the perturbation potential caused by a thin 
wing in a uniform flow of Mach Number J/). The 
orientation of the coordinate axes relative to the free- 
stream direction and the wing is shown in Fig. 1. 

Supplementary requirements on ¢ require it to vanish 
on the Mach envelope from the wing leading edges, 
and to satisfy the boundary condition on the wing, S,,, 
which requires that the flow downwash angle be coinci- 
dent with the wing angle of attack distribution a(x, y); 
that is, 


—¢.(x, y, 0) = a(x, y) on Sy, (2) 
Then, if s = (x, y) is the wing shape, 
—a(x, y) = 2,(x, y) 

or —a(x. y) = ff a(x, + g(y) on S, (3) 


Here the integral is an indefinite integral, and g is an 
arbitrary function. 

The wing problem formulated above has been solved 
for a general plan form (see references 13, 14, and 15, 
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for example), and for our purpose we can consider 
the wing pressure coefficient to be known as a functjoy 
of the given wing angle of attack—that is, we have 


Crp = C, [a(x, y)] (5 


The coefficients of drag, lift, and pitching moment 0 
the wing of area S,,, and mean aerodynamic chord C ax 
given by 


SuCp = a(x, y) dx dy 
Su 
S,Cr = 2 C, dx dy (6 
Su 
2 f { xC, dx dy 


In the expression for Cy, the origin of the coordinate 
is at the leading edge of the root chord. 

The problem is now to determine an angle of attack 
distribution @ = app. (x, y) for a wing of given plan 


form and flight Mach Number which yields a minimum | 


value of Cp for a prescribed value of C, = C,, and 
Cy = Cy, and possibly additional geometrical con- 
straints. 

For the solution of this problem we first assume that 
the a may be represented by the series 


a(x, y) = > Cyx* (7 
i 


where the coefficients C;; are now to be determined. 
The requirement that the wing yields the prescribed 
values C, and Cy is next expressed by the conditions 


g(a) = C, — Ci(a) = 
g(a) = Cy — Cuy(a) = Of 


We then introduce two new parameters \; and do, the 
so-called Lagrangian multipliers, and form the function 


F > Cp + Aigi(a@) a) 


where a = a(C;;) from Eq. (7). The extremum of Cp 

is obtained when the function F fulfills the system oi 

equations 

OF/O\, = g = 0, & = 1,29 

Note that the last two equations are just the constraints 

given by Eq. (8). Eqs. (9) now form a determined 


system of linear equations for the determination of the | 
C;; and the d; and \, when a partial sum is taken in | 


Eq. (7). 


Any geometrical constraints on the wing shape will | 


decrease the degree of freedom by which the wing shape 
may be permitted to vary for the drag minimization. 
This is effected by the additional relationships between 
the coefficients C;; required by the geometrical con- 
straints through Eg. (3). These relations must be in- 


serted into Eq. (9) to decrease the number of available | 


coefficients C,;. Note in Eq. (3) that the arbitrary 
function g may be used to shear the wing to obtain, for 
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OPTIMUM THIN 
instance, one straight control surface hinge line on the 
wing without constraining the C,,’s. 

The minimum value of Cp can then be computed 
using the value of the Lagrangian multipliers by the 
simple expression (see next section) 


NCL + (10) 


—2Cp.., = 


Aside from the difficult mathematical question of 
completeness of the set of basic angle of attack distribu- 
tions, as mentioned in the Introduction, there is addi- 
tionally the practical question of the “‘naturalness”’ of 
the chosen set of a@ functions as manifested by the 
rapidity of convergence of the process. In general an 
indication of the convergence is obtained by taking 
more and more terms in the a-distribution, and noting 
the rate at which the sequence of resultant drag values 
approaches a limit. Here there is, of course, still the 
possibility that one does not obtain the smallest limit, 
since the lack of completeness of the basic angle of at- 
tack distribution in Eq. (7) may act as a constraint on 
the generality of the a@ distribution. To check this 
possibility, the minimum values of the drag are also 
computed by an independent method to be described 


in the following section. 
HARMONIC METHOD 


Here we shall consider the harmonic method derived 
independently by Nikolsky’ and Kogan‘ using a control 
surface approach and by Graham’ using the ‘“‘combined 
flow’ concept. These methods give the minimum value 
of the drag independent of the set of basic angle of 
attack distributions required in the method of the pre- 
vious section. However, this procedure is somewhat 
inconvenient to use in finding the necessary wing shape 
to produce the minimum drag since either its domain of 
applicability is confined to the rear Mach envelope or 
it involves directly only the combined flow field. Thus 
the potential theory approach and the direct approach 
of the previous section supplement each other. 

We shall first give a brief description of the ‘‘com- 
bined flow” approach to the potential theory method 
and then give a derivation using the control surface. 

The basis of the combined flow approach is the opti- 
mum drag criterion derived by Jones" pertaining to a 
wing which yields minimum drag due to lift (vortex 
plus wave drag) when the lift, wing plan form, and 
flight Mach Number are prescribed. To illustrate the 
result, consider a symmetric diamond plan form with 
subsonic edges. (Jones’ result is applicable, of course, 
to any plan form and to both subsonic and supersonic 
Mach Numbers.) We now define the 
lorward flow when the free stream approaches the wing 
irom the left. 
lorward flow, we will obtain a given pressure distribu- 
tion. We now keep the position of the wing fixed and 
reverse the free-stream direction such that it now ap- 


free-stream 


For a given angle of attack a, of the 


proaches the wing from the right; we then determine 
the angle of attack a, of the reverse flow which will 
vield the same pressure distribution on the wing as in 
Jones found that, for minimum drag, 


the forward flow. 


LIFTING 
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Domaia of the harmonic problem 


Fic. 2 


it was both necessary and sufficient that the sum of 
the angle of attacks in the forward and reverse flows 
be a constant—1.e., ay + ag = constant, on the wing 
plan form. 

Jones’ criterion for minimum drag in the case of 
wings with subsonic edges did not include the con- 
sideration of the leading-edge suction force in the ac- 
counting of the drag.** Strictly speaking, in any vari- 
ational procedure using Jones’ criterion, admissible 
angle of attack variations must be restricted to those 
which yield finite pressures at the wing edges. 

Furthermore, Jones found that the optimum drag 
in either the forward or reverse directions could be ex- 
pressed in terms of quantities in the combined flow 
field. 

Graham, using the combined flow concept, intro- 
duced physically artificial forward and reverse flows 
composed of line sources and vortices which gave the 
required combined flow field for minimum drag. Here 
the requirements for the optimum combined flow field 
include the fulfillment of Jones’ criterion on the wing 
surface, and the maintenance of the proper zone of 
silence and the proper vorticity distribution in the 
wakes far upstream and downstream of the wing. The 
resultant combined flow field is then two-dimensional 
within the domain formed by the envelope of the Mach 
cones from the wing leading edges in forward and re- 
verse flows. The combined flow in this region can 
thus be formulated as a two-dimensional potential 
problem from which the optimum drag can be deter- 
mined. 

Jones’ criterion, therefore, is basic in the procedure 
of Graham. The optimum drag derived by Graham 
then is correct only where the wing pressures are 
bounded and leading-edge suction is absent. 

The method of Nikolsky and Kogan uses a control 
surface formed of the forward and reverse Mach en- 
velopes, as in Fig. 2, and computes the aerodynamic 
forces and moments in the forward flow by using the 
usual conservation laws. The variational procedure 
with Lagrangian multipliers of the previous section is 
next applied using these expressions for the drag, lift, 
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and moment; one then obtains a potential problem 
identical to that found by Graham, but in terms of the 
flow variables on the rear Mach envelope. Here, of 
course, the potential of the forward flow is identical to 
that of the combined flow since the reverse flow poten- 
tial is identically zero there. 

The use of the control surface has the advantage over 
the combined flow approach in that any leading-edge 
suction force is automatically included in the drag ex- 
pression. * 

In the case where only the lift is prescribed we see 
that the identical harmonic problem is_ obtained, 
whether it is derived using the combined flow approach 
where the leading-edge suction is assumed to be absent 
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or by the control surface method where the leading. 
edge suction is automatically included. Thus, it may 
be concluded that the optimum wing shape in the cas } 
of subsonic leading edges must be such that the pressures 
are bounded at the leading edges and that no leading 
edge suction arises. Moreover, the Kutta condition 
must be fulfilled at subsonic trailing edges. 

In the following we shall rederive the harmon, 
method based on the control surface approach wher 
both the lift and pitching moment are prescribed 
The general orientation of the wing and the controj 
surface is shown in Fig. 2. The conservation laws 0; 
momentum and moment of momentum yield'® th 
following equations for the coefficients of drag, lift,"and 
moment: 


Cp = (1/0 (v? + w? + u} — 2v{cos (», y) ‘cos x)] — 2w[cos (», 2) cos (», x) ]}) dS 


Cy = Seti f }w — u [cos (r,s) cos (vy, x) |! dS 
s 


Cu = (2/CUS,) + v [cos y) ‘cos (v, x)] + w [eos s) ‘cos (vy, x) + 


Here 
(u, v, w) = components of perturbation velocity 
vector 
= — 1)'/*, My being the free-stream 
Mach Number 
Uo = free-stream velocity 
cos (v, x), 


direction cosines of the interior normal 
cos (vy, = 
of the rear Mach envelope 


cos (y, 3) 

dS = projection of an element of surface 
area of the Mach envelope on the 
plane transverse to the flow direc- 
tion 

S = projection of rear Mach envelope on 
the same transverse plane 

Sis = wing area 

C = wing mean aerodynamic chord 


In Eq. (11) the surface integral is to be taken over 
only the rear Mach envelope since the perturbation 
velocities vanish just upstream of the forward Mach 
envelope. 

Next let x = f(y, z) be the equation of the rear Mach 


envelope. We have then the following relations: 


cos (vy, x)/—1 = cos (», y) f, = cos (», 


We next introduce into Eq. (11) the change of notation 


(12) 


&(y, z) = ¢lf(y, 2), y, 2 (13) 


* On the other hand, for wings with subsonic trailing edges, 
as pointed out by Ursell and Ward," the occurrence of infinite 
velocities will locally invalidate the linear theory at the side edges 
of the trailing vortex sheet where it intersects the rear control sur- 
face. 


x} —w + u [cos 2) cos (», x)]!) ds | 


with = + ¢. + ¢,/: 


We obtain now the simple expressions 


l 
Cp (,° ®.*) dS 


C | US,» dS (14 


9 


CUS. + 2(,f, + 


Using these expressions we again form the function 
f= + igi + Avge 


as in the previous section, and then set the first vari- 
ation of Fequal to zero. The resulting Euler equation” 


and the natural boundary condition are 


® = ).(2f, + 2V°f) in Sr] 
and ond (15 
=r, + Af 


where the domain Sp and the contour C,; and (>, are 
shown in Fig. 2. 

The function ® is unsymmetric with respect to the: | 
coordinate. Eq. (15) is now the required boundary 
value problem to be considered. We shall carry out 
a further transformation to obtain the form of the 
problem obtained by Germain—that is, we put 


P(y, 2) = Ale + xily, + + (16 


The X,; and X, are defined by the boundary value prob- 
lems 
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Fic. 38. Plan-form shapes (sonic leading edge 


= O in Sp, 
Xi = onC, (17) 
X12 = 0 on 

and = O in 
X = -sf ond, (1S) 
= on cA 


The boundary value problems given by Eqs. (17) and 
(18) can be solved by standard methods (see the follow- 
What remains is the determination of the 
These can be evalu- 


ing section). 
Lagrangian multipliers and A». 
ated by inserting the resulting value of ® into the ex- 
pressions for lift and moment and inserting these, in 
The latter 
equations then give two linear simultaneous equations 


turn, in the constraint equations (S). 


for and Xe. 

A simple well-known expression tor CDent. given pre- 
viously in Eq. (10), can be shown following Heaslet.' 
The procedure is to transform the various surface in- 
tegrals in Eq. (14) for the lift, drag, and pitching mo- 
ment by means of Green's formula involving the La- 
placian and then, using Eq. (15), to identify the various 
terms in the drag expression with those for the lift and 
pitching moment. 


SOLUTION OF THE HARMONIC PROBLEM 


For specific plan forms there are special procedures 
by which simple analytical solutions are possible (see, 
for example, references 9, 12, 20, and 21). For the case 
of a general plan form, however, we shall use the follow- 
ing method. 

Since the procedure for calculating yx; and x» is the 


same, we shall consider only the latter. We first put 


= x" x” (19) 


where the harmonic functions x‘!’ and x‘? fulfill the 
boundary conditions 


qd) — — of (2) — 0 
x f, x on Cy (20) 
= = —x, on 
For the determination of x“) we set 
x) = > A,P,(y, 2) (21) 
n 


Here 


P,, = esch n[cosh n(y ‘y,,) sin 2 (2 Zn) + 


Where y,, and z,, are the maximum values of the y and 


L 
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zs coordinates on C;. In Eq. (21) the coefficients A, 


are now determined such that the condition y‘') = —2/ 
is fulfilled on C;. The form of P, 
trial and error process of evolving a well-behaved ma- 
by the 


was evolved by a 


trix in the determination of the coefficients A, 
collocation method. Additionally P,, must be symmet- 
ric with respect to y and antisymmetric with respect 
to 2. 

For the determination of y‘*’ we shall first form a set 
which fulfills 


These particular solu- 


of harmonic functions Q;(y, 2), each of 
the condition Q; = 0 on C). 
tions will then be superimposed such that the remain- 
ing condition = —x.'' on Cy is fulfilled. 

For the function (Q we set 
vy) B, 


n 


Oy, 2) = 3; P, (22 


where the function I'(y, 2; y,;) represents a pair of finite 
vortices of equal and opposite strength situated on the 
wing slit C. symmetrically about the z axis at points 
(y,, 0) and (—y,, 0). The family of solutions Q, is then 
formed by placing the vortex pair at different points 


y = +y, along the slit C.. The coefficients B,'" are 
now determined from the condition Q; = 0 on C. 
Lastly, we have 

(23 


x2? = CO 


where the coeflicients C, result from the fulfillment of 
the remaining condition on the normal derivative along 
C.. The sum of Eq. (23) and Eq. (21) then will give 
our desired solution. 


RESULTS OF THE CALCULATIONS 


In order to illustrate the direct method, a delta wing 
and a sweptback wing with sonic leading edges (see 
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Fig. 3) were computed with the constraints prescribing 
either the lift alone, lift and pitching moment, or lift, 
pitching moment, and two straight hinge lines. 
Ten terms of the angle of attack distribution of Eq. (7) 
were used for the calculation, and the static margin 
in all cases was taken as 14 per cent. Thus, the trim 
drag is zero at the design lift coefficient and Mach 
Number. The resulting wing warp for the various 
cases is shown in Fig. 4. It is to be remembered here 
that the various chordwise shape distributions may be 
displaced with respect to each other within reason- 
able limits according to Eq. (3). The corresponding 
values of the drag due to lift are as follows: 

= (Cp/BCz?) 


Sonic Sonic 
Constraints delta wing swept wing 
Lift 0.2253 0.1835 
Lift, pitching moment (static mar- 
gin: 14°%) 0.2378 0. 1862 
Lift, pitching moment, two straight 0.2945 0.1965 
hinge lines (static margin: 14°,) 
Flat (plane) wing 0.250 0.234 


In order to check the sufficiency of the ten term a 
distribution, the same plan forms were also computed 
by the harmonic method for the case in which only the 
liftis prescribed. The results for the sonic edge case, as 


well as those for the subsonic and supersonic edge cases 
are shown in Fig. 5. The boundary conditions were 


satisfied at 20 or more spanwise points. This gaye 


answers that were within 0.5 per cent of the cases re. 
ported in references 9, 12, and 27. 

A comparison between the results of the harmonj 
and the ten term direct methods show excellent agree. 
ment. Some results obtained by Doris Cohen? for g 
delta plan form with sonic and subsonic edges are 
plotted in Fig. 5(a). Here six quasi-conical terms jy 
the a-distribution were used. The results of Doris 
Cohen include leading-edge suction and are in good 
agreement with the harmonic method (except at very 


low 8.R where a slight discrepancy exists) despite the \ 


fact that the latter theory indicates the absence oj 
leading-edge suction force for the optimum case. The 
quasi-conical a@ distributions in this case would not be 
the natural ones because of the infinite pressures arising 
for each of the basic distributions. 

In reference 28 Germain and Gibault have plotted the 
drag reduction possible for a delta wing with subsonic 
and supersonic leading edges. A discrepancy between 
their results and ours is apparent at 8.R < 2.5. We 
have carefully checked our inputs to the digital pro- 
gram without finding any errors. An increase in the 
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Fic. 5(a). Drag due to lift—delta wing. 


C, CONSTRAINED 


| | | 
| | | | 
FLAT PLANE) WING | 
OPTIMUM WARPED 
WING | 
NOLE | 
SUCTION ~ 
28 T | 
| es 
FULL Le | | 
SUCTION 
| 
16 | 
o 2 3 4 


Fic. 5(b). Drag due to lift--sweptback wing. 
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OPTIMUM THIN 
number of terms used from 20 to 30 shows only a very 
slight decrease in the value of Cp 8C,° reported in 
Fig. 5(a). Moreover, the drag values for 15, 20, 25, 
and 30 term calculations indicate convergence of the 
procedure. 

In Figs. 6 and 7 we show the results for additional 
plan forms using the harmonic method where only the 
lift is prescribed. For supersonic edge cases the results 
generally indicate that the drag reductions are in- 
creased with wing sweepback and with taper ratio. 
Furthermore, the drag reduction for a given plan form 
decreases as would be expected for increasing Mach 
Number. 

Lastly, in Fig. 8S we show the possible gain in the lift- 
drag ratio due to wing warp for a typical sweptback 
wing. Here the maximum lift-drag ratio is given by 


[L Pines (1 2) [Cp, (Cp 


where Cp, is the zero lift drag coefficient excluding drag 
due to camber and twist. Here it is seen that the in- 


crease in lift-drag ratio due to warp may be substantial. 


CONCLUDING REMARKS 


Supersonic interference effects furnish a fertile field 
to seek improvement in the performance of high-speed 
aircraft and missiles. In the present paper the simplest 
of these effects, namely, the supersonic interference 
between lifting elements of a wing surface has been 
investigated. In a practical configuration, however, 
one must additionally consider effects arising from 
wing-body interactions. 

With the harmonic method it was shown that, for 
the optimum wing with subsonic edges, the pressures 
were bounded, and no leading-edge suction force arises. 
Practically speaking, the absence of the leading-edge 
suction force has an additional beneficial effect by 
directing the wing leading edge into the oncoming 
stream, thereby minimizing any leading-edge boundary- 
layer separation effects. 

The ten term quasi-conical representation of the a 
in the direct method was found to be adequate to deter- 
mine the minimum drag as well as the required wing 
warp for the plan forms considered with the sonic lead- 
ing edges. 

For the case of wings with subsonic edges, however, 
the quasi-conical basic angle of attack distributions 
would not be the natural ones. The basic a-distribu- 
tions used, for example, by Tucker,”? where only finite 
pressures resulted, would be more appropriate for the 
treatment of wings with subsonic edges. 

Lastly, caution must be exercised in using the opti- 
mum drag values from the harmonic method, since 
unrealistic wing shapes may be required to attain these 
values for some plan-form shapes. 
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Angle of Attack Convergence of a Spinning 


Mis Sile Descending Through the Atmosphere | 


HERMAN I. LEON* 
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SUMMARY 


In many missile designs, aerodynamic damping is negligible at 
altitudes above that for maximum heating. In addition, vari- 


ations in velocity and flight path angle are small. Utilizing these 


approximations, the angular motion of a spinning missile de- 
scending through an exponential atmosphere is derived for small 
angles of attack. 

A certain amount of angle of attack convergence is experi- 
enced by a nonspinning missile, attributable to the exponentially 
increasing air density. It is found that for a symmetrical missile 
spinning slowly, and having no initial angular rates, the envelope 
of angle of attack near maximum heating is increased over that 
for nonspinning re-entry by the ratio 


[( tanh (ak;/2)|!/? 


where &; is a spin parameter proportional to the spin rate. 
For a rapidly spinning missile under the same conditions, the 
ratio becomes 
(a4 2)! ks?) + tk. 
where k» is a static stability parameter proportional to the mo- 


ment coefficient. 


SYMBOLS 


Cp = aerodynamic drag force coefficient 

= aerodynamic lift force coefficient 

Cy = aerodynamic moment coefficient 

Cap = (0C4/0B)p =0 

g = gravitational acceleration 

Ty, Jy, /z = moments of inertia about the pitch, yaw, and 
roll axes, respectively 

I = Ty or ly 

= -1 

- = Bessel function of the first kind 

= sin? yx) = static stability 
parameter 

= Vr sin ye) = spin parameter 

l = reference length 

m = mass 

ae Pe = components of the angular velocity along the 
spin, pitch, and yaw axes, respectively 

R = radius of the earth 

Ss = reference area 

J = velocity 

XY, ¥,2Z = missile pitch, yaw, and spin axes, respectively 

VY’, 3’, Z’ = inertial axes having Z’ aligned with the flight 
path 

= real variables 

y = altitude 

a = complex angle of attack = 6e~** 

8 = air density exponent factor 

SY = flight path angle from local horizontal 

v = complex number 
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air density 


i 

v,0 = Eulerian angles relating Y, Y, Z and Y’, 
Z' (see Fig. 1) 

= -—1)/I 

w = complex tumble rate = g + ir 


Subscripts and Superscripts 


0 = reference, constant 
E = re-entry 
= time derivative 


INTRODUCTION 


ARE many considerations in the re-entry 
design of a missile that require prediction of its 
angular motion as it descends through the atmosphere. 
Problems relating to stability and to the effect of angle 
of attack on loading and heating in missiles during re- 
entry have been given considerable attention.’ *° 
Previous work, however, has been confined to non- 
spinning bodies entering the atmosphere in a yawed 
or pitched attitude. A spinning missile could be ex- 
pected to exhibit a different dynamic behavior. This 
paper extends the angle of attack analysis to the case 
of a re-entry body entering the atmosphere with an 
angle of attack and spinning. 

For present missile systems in which re-entry is a seri- 
ous problem, the required angle of attack convergence 
is specified almost completely by heating and loading 
effects. It has been shown? that maximum average 
heating occurs when the velocity has decreased to 
0.72 Ve, and maximum stagnation point heating at 
0.851. Also, maximum axial loading occurs at 0.61 
Vy, while maximum normal loading occurs at this ve- 
locity, or at a higher velocity if aerodynamic damping 
is appreciable. The important region for angle of 
attack convergence, then, is the altitude range above 
that corresponding to these velocities. 

The velocity change during re-entry is derived in the 
Appendix. The form of the velocity equation is such 
that the velocity change is rather abrupt with altitude, 
changing from 0.95 V, to 0.61 Ver in about 50,000 
ft. It isin this 50,000-ft. altitude range that the angle 
of attack must be confined to a relatively small value. 
The convergence, therefore, must occur at higher alt- 
tudes that is, before the velocity has changed, say, 
5 per cent from re-entry velocity. Therefore, a solt- 
tion to the re-entry problem that postulates constant 
velocity will be quite meaningful. The accuracy of 
such a solution will probably be at least of the same 
order as that of the aerodynamic coefficients in the 
very high altitude ranges of interest. In this analysis, 
a constant velocity is postulated. 
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Z' FLIGHT PATH DIRECTION 
Z NOSE CONE SPIN AXIS 


@ MAGNITUDE OF ANGLE 
OF ATTACK 


Fic. 1. Definition of Eulerian angles 6, ¢, ¥ relating the n 
cone axes Y, Y, Z to the reference coordinate system Y’, }’, Z 


aligned with the flight path. 


se 


In the altitude ranges of small velocity change, the 
flight path angle is practically constant (see Appendix), 
except for grazing re-entry of a few degrees or less with 
the local horizontal. A constant flight path angle is 
postulated in this analysis. 

The effect of aerodynamic damping has been analyzed 
for nonspinning re-entry bodies.*° These results show 
that convergence is improved by an amount that de- 
pends on the aerodynamic coefficients and further, 
that this damping effect is a rather abrupt function of 
altitude. Terms in the differential equation that arise 
if the velocity is allowed to vary also contribute to this 
damping term. 

The form of the solution (see Appendix) is such that 
the damping effect is negligible down to altitudes cor- 
responding to 0.95 Ve for blunt bodies. For long, 
slender bodies, the damping effect may be considerable 
at thisaltitude. In this analysis, aerodynamic damping 
and the comparable term arising from differentiating 
the velocity are neglected. Therefore, the solution pre- 
sented here is accurate for blunt bodies, and provides a 
conservative limit on angle-of-attack convergence for 
the long slender shapes. 

The variation of air density will be considered an 
exponential function of altitude, 


poe”? (1) 


where p, and 8 are constants. 

Because of the simplifying assumptions used in the 
analysis, no effect of spin on stability can be revealed 
here 


EQUATIONS OF MOTION 


Consider a missile descending through the atmosphere 
and having a center of gravity trajectory of constant 
velocity and constant path-angle. Let its angular 


481 


orientation be described by a set of Eulerian angles as 
defined in Fig. 1. 

For a symmetrical missile having no products of 
inertia, /y and /y are equal. (The single symbol / will 
be used.) In the altitude region of interest the only 
moment is the aerodynamic restoring moment, which 
is proportional to the angle of attack for small angles. 
The rotational motion is described by Euler's dynamic 
equations 


0 = 


(1 2)pV*S/Cy, 6 cos Y = Ig — UI — Iz)pr (2) 


—(1 2)pV*S/Cy, 6 sin Y = + UI — Iz) pq 

The first equation in (2) integrates immediately, and 
The 
remaining pair of equations can be combined in terms 


shows that the roll rate must remain constant. 


of a complex angular rate (tumble rate) as 


= 


(3) 
Eq. (2) becomes 


p( 27) Cu, = — p= (A) 


where \ = (Jz — 1) J. The relationship between the 
body rates and the Eulerian angle rates is given by* 


p=¢cos6+y | 
q = @sin@sin + 6 cos (5) 
r = @sin@cos — 6sin 
From Eq. (5) we have 
w= (6+ i(p, — tan dle” (6) 
Let us define a complex angle of attack as follows: 
a =e" (7) 


This angle of attack is related to the complex rate (for 
small a) by 


a+ Ip,a (S) 


Rewriting Eq. (4) in terms of this new variable, we have 
@+ i(1 — A)pa + 


= 0 (9) 


For an exponential atmosphere, and in terms of alti- 
tude, 


(d?a + i8ks|(A — 1) (A + 1)] (da/dy) + 


B? | + + 1)?]f = (10) 
where dy = —Vesin ye dt 
The constant, 

ko = sin? ye (11) 


termed the static parameter, determines how stable 


the missile will be in the nose-first attitude. The 
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| 
constant, a(y) = V xk; sinh rk; exp | —(i/2)k3[(A — 1) + 
ks = plz sin ye (12) — 1)] By} + x 


will be called the spin parameter. 
The solution to Eq. (10) is given by 
a = IBY [4 (2V 4 
BJ~u, (2W (13) 
The coefficients A and B are determined by the initial 
conditions of the missile. 
The Bessel functions of imaginary order are defined 
by the usual Bessel series expansion' 
x 27 +ik 
(14) 


x )j 
Jax) = ; 
+1 + ik) 
o.* = 0; and 


At the initial condition point of y = 
lim = [(x 2)"/T(1 + ik)] = 

—i(x/2)* V (sinh wk) mk (15) 
The general solution of Eq. (13) evaluated at infinite 


altitude becomes 


(sinh wk;)/ wk; X 


lim a (y) = 


2 Bp, — ike 2 (16) 
In terms of time, 

lim a(t) = —iV (sinh wk) wk; X 


The constants A and B can be determined by solving 
Eq. (4) and (8) for the case of p = 0 (no atmosphere). 
This solution, which is for infinite altitude, is 


a = age‘? + + A)] X 


(18) 


where ag = 6» is the initial angle of attack and we = 


de + ire is the initial tumble rate. The coefficients 
become 

A = [wr /p,(1 + d)] V rk? /sinh wk; 


Works sinh ahs 


The complete solution is 


j=l 


we mie 


For an imaginary order, 


J x(x) ~ V2 xx }1 — [(4k? + 1)/16x?]! cos }x — (w/4) — tan—! [(422 + 1) Sx] — 


The first significant terms in x have been retained in order to estimate the accuracy of the approximation 
nitude of the cosine of a complex angle is bounded as follows: 


— 1°] — — (4 — 


Jing (2W + + + d)]} x 


NATURE OF THE SOLUTION 


The initial conditions, Eq. (18), can be interpreted 
on a physical basis. 
at the same velocity but ahead of the missile, and ob. 
serving the spin-up procedure, considered here for con- 


Consider an observer traveling 


venience as an instantaneous change KC.2 by the use 
of solid fuel rockets of short burning time). Other 
methods of spinning produce an equivalent motion 
after the spin torque has been exhausted. In his non- 
rotating coordinate system, he would observe the 
angle of attack behavior after spin-up as 


ac?" = ap +i [we pl + — (18a) 


The tip of the missile would exhibit a circular motion, of 
diameter 2we p,(1 + A), and frequency p,rA 27 eps, 
passing through the angular position of spin-up, az, 
at some point on the circumference. The circle diam- 
eter would be at right angles to the direction of tumble. 
the would lose its initial 
significance and become a circle of the same diameter 
and frequency, but it would be displaced from the zero 


To this observer, motion 


angle of attack position by an amount 
ag + [twe p(1 + 


This total average angle of attack gives rise to the 
/J_, term in the general solution of Eq. (20). The 
circular motion attributable to initial tumble rates 
gives rise to a nutation-like motion described by the 
J i, term of Eq. (20). 


(18b) 


RELATIVE CONVERGENCE 


Low Spin Rate 


The primary interest here concerns the amount of 
angle of attack convergence obtained during the 
descent through the atmosphere down to altitudes 
near maximum heating. For many missile con- 
figurations, the argument of the Bessel function is ex- 
tremely large when evaluated for these relatively low 
altitudes. Therefore, the asymptotic expansion for a 
Bessel function of large argument is appropriate. The 
expansion for small values of pv is given by’ 
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(2j) 126/02 


(—) [4v? — 1°] [40 — 3]... [407 Gj (21) 
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ATTACK CONVERGENCE 
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A SPINNING 


OF 


(23) 


< cosh 


The envelope of the complex angle of attack for large values of the Bessel function argument and low spin rates be- 


comes, neglecting the lower order terms in x, 


Env fa(v)} ~ V rk; sinh zk; * V eV bo) + A)] cos [E(y) — tks 2)] + 


where £(y) is real. For the case of a missile spinning 
slowly without wobble (we = 0), and entering the 
atmosphere with an initial attitude error, ag, the angle 
of attack envelope reduces to a comparatively simple 
expression, 

tk:; 2 


tanh rk; 2 


Ape 
V be 
In this case Eq. (25) reduces to the planar solution under 
the same restrictions when k; and p, are set to zero. 
The point of interest is the relative convergence of a 
slowly spinning missile without wobble over that of a 
The in- 


Env la(y)} (25) 


nonspinning one under the same condition. 
crease in the envelope of the angle of attack attributable 
to low values of spin is 


tanh rk; 


Env a(y)} 


(26) 


Env 

The most significant term neglected (compared to 
unity) in this derivation is 

(4k;? + 1) (27) 

Neglecting this term is automatically justified by the 

more stringent condition for the expansion in Eq. (21); 


x exp (ik Vi¢x k?) exp (— Wx) 


2 
(QE 


+ [iwe p,(1 + cos [E(y) + mar (24) 


4 ‘ f 1,2 
viz., ky 2(kee™ 


of spin is considered below. 

Eq. (26) is plotted in Fig. 2... For spin value param- 
eter k; > 1.5, the hyperbolic tangent term is nearly 
unity, and the relative convergence depends on the 
The relationship shown 


Convergence for larger values 


square root of the spin rate. 
in Fig. 2 may be considered as the penalty factor on the 
angle of attack envelope for the use of slow spin for the 
zero wobble re-entry case. 

The case of an initial tumble rate is inconvenient for 
comparison because of the infinite altitude boundary. 
In the planar case (p, = O) an initial tumble would 
generate a random angle of attack at any finite alti- 
tude, and hence would violate the small-angle condition 
for this solution. A finite altitude boundary may be 
used, but the solutions become considerably more com- 


plicated. 


High Spin Rate 


When the order of the Bessel function is large then 
the asymptotic expansion for large argument is con- 
siderably different. For the case of an imaginary order 
and a real argument, both of large magnitude, the 
Meissel expansion is suitable :* 


(28) 


J x(x) 
e"T (1 + + x? V1 + x? 
where 
W, = l | 2 — 3x?/k? 2] 4 xi kt — 4x? R? l 
+ 1512 x2/k? — 3654x1/k* + 375 | 32 x? — 288xi/kt? + 232 — 13x5/R5 

9 9\9/2 16 + + 
(1 + 128 kt (1 + x? 


Inspection of the series for WW, shows that for large 
values of k and x, the real part of IV’; is negligibly small 
for both large and small values of the ratio x k. The 
magnitude of the asymptotic expansion becomes 


sinh l ) 
(a) ( wk + x?/k? 


The asymptotic envelope of the complex angle of at- 
tack for large spin rates can be evaluated from Eq. 
(20). For the case of a missile spinning rapidly with- 
out wobble (we = 0), and entering the atmosphere 
with an initial attitude error, az, the angle of attack 
envelope reduces to 


QE 


Env } a(y)} 


The relative convergence of a rapidly spinning missile 
without wobble with respect to a nonspinning missile 


1/2 1 4 


CONCLUSIONS 


becomes 


Env } a(y)} 


Env a(y)} 


The angle of attack envelope is increased during re- 
entry for a spinning missile over that for a nonspinning 
one. 

For the case of an initial attitude error and no wobble 
or tumble rate, the angle of attack envelope is increased 
over nonspinning re-entry by an amount dependent on 


. 
| | 
= 
) 
q 
| 
. . 
| 
= 
(22) 
mag- = 


484 JOURNAL OF 
— | | 
\22 
Env tanh ks | | 


a =complex angle of attack 
4 4 


k, = spin parameter 


Z 


Relative Convergence 


' | 
Oo ! 3 4 


Spin Parameter kz 


Fic. 2. Relative convergence of the angle of attack envelope for 


spinning re-entry (low spin case). 


the spin rate, altitude, and trajectory and missile 
parameters. 

The effect of initial tumble rates at spin-up generates 
a wobble, which gives rise to a nutation during re-entry. 


APPENDIX 


The trajectory equations for the center of gravity 
of the missile are given by” 

V = —(CpS’m) p(V?/2) — g sin 
(32) 

y= {[V2/(R + y)] — gf y)/V]\ 
The second equation can be used to predict the path 
Note that the 1*/(R + y) term is 
A conservative 


angle change. 
equal to g only for satellite velocity. 
estimate of y is simply to neglect 1° (Xk + y) with 
respect to g. Therefore, change in path angle in the 
altitude range where velocity change is small is approxi- 
mately 


dy/dy & tan yz) (33) 


For missiles that are re-entering the atmosphere at 
other than grazing incidence, the flight path angle 
under the above conditions is less than a degree. 

The first equation of (32) can now be sclved for the 


case of constant path angle. Neglecting the gravity 


Change 
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term* and for an exponential atmosphere and constant 
Cp, the velocity is found to be! 


= Ve exp | —(Cpp,S/28m sin yz) (34 
The dynamic pressure becomes 


qg = (1/2) pV? = (1/2)p,V exp —B8y — 
(Cp 3m) (35 


The altitude for maximum qg occurs wheu the exponent 
in Eq. (34) is —(1/2). The altitude at which the 
velocity is, say, 95 per cent of the re-entry velocity js 
found to be about 50,000 ft.t above the altitude of maxi- 
mum g. Therefore, the assumption of constant ve- 
locity is reasonable in the altitude range starting about 
50,000 ft. above the altitude of maximum gq. 

Aerodynamic damping gives rise to a factor in the 
planar angle of attack solution *° 


exp | —(p,4,/48m sin yz) X 
[Cp = + (Cy, + Cy,)ml* 1] e (36 


At maximum gq, the term becomes 


jl a Cu, + Cue 
exp = (54 
\4 Cop D I 


At 50,000 ft. above maximum g, this exponent is re- 
duced by a factor of 10. At this altitude, the damping 
term will contribute a change of, say, 5 per cent or less 
when the magnitude of the aerodynamic term in 
brackets in Eq. (37) is greater than 2. This analysis is 
restricted to such missiles. 
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SUMMARY 


\ survey of existing analytical treatments of the supersonic or 
hypersonic blunt-body problem indicates that none is adequate 
for predicting the details of the flow field. Reasons are given for 
the failure of various plausible approximations. A numerical 
method, which is simpler than others proposed, is set forth for 
solving the full inviscid equations using a medium-sized electronic 
computer. Results are shown from a number of solutions for 
bodies that support detached shock waves described by conic 


sections 


SYMBOLS 


B = bluntness of conic section |see Eq. (1 
= 1—B 
G = pressure coeflicient referred to free-stream condi- 
tions 
pip 
WV = free-stream Mach Number 
D = pressure referred to po Vx? 
R = nose radius 
= coordinate normal to free stream 
= components of velocity in &, directions 
\" = velocity vector referred to free-stream speed 
\ = coordinate in free-stream direction 
= adiabatic exponent 
6 = ] + +1) 
A = standoff distance of shock from body nose 
gn = orthogonal curvilinear coordinates |see Eq. (2 
Ag, An = mesh widths in numerical computation 
y = 0 for plane flow, 1 for axisymmetric flow 
p = density referred to free-stream value 
Vv = stream function 
w = modified stream function {see Eq. (11 
6 = polar angle 


= value on shock wave 
b = value on body 
= free-stream value 
= value at stagnation point 


INTRODUCTION 


P" TOGRAPHS of the detached shock wave that forms 
ahead of a blunt body in supersonic flight show a 
beautiful simplicity that has in the last decade lured a 
number of aerodynamicists into trying to predict its 
shape. Recently there have arisen more urgent rea- 
sons for finding the complete flow field (particularly in 
hypersonic flow), and a considerable fraction of the 
theoretical aerodynamicists in the United States, 
England, and (apparently) Russia have become in- 
volved with the problem. The attacks have been 
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mainly analytical, though within the last year elec- 
tronic computers have been brought to bear. 

The present paper aims, first, to evaluate the mass 
of existing analytical treatments and, second, to put 
forth a new and relatively simple numerical procedure. 
The survey will show that existing analytical methods 
are generally inadequate for predicting the details of 
flow near a blunt nose. With the numerical method, 
accurate solutions have been carried out for members 
of a one-parameter family of plane and axisymmetric 


bodies. 


(1) REVIEW OF EXISTING ANALYTICAL SOLUTIONS 


Existing analytical attacks on the supersonic blunt- 
body problem fall mainly into four categories: (1) 
Potential flow approximations.'~* (2) Taylor series 
11 » 


expansions from shock.* Incompressible ap- 


proximations.'*~" (4) Newtonian approximation and 


improvements thereon. 
The only other noteworthy approach is that of ex- 
which gives the 


plosion and similarity solutions*® 
asymptotic flow field downstream on a slender blunt- 
nosed body. This elegant and useful theory is not in- 
tended to treat the vicinity of the nose and will, accord- 
ingly, not be considered here. The four methods listed 
above will now be discussed in succession, drawing 
where necessary on the accurate numerical solutions 
described in section (2) as a standard of comparison. 
Other will be 


section (3). 


numerical methods*~** discussed in 


(1.1) Potential Flow Approximations 


The oldest attack on the blunt-body problem con- 
sists in attempting somehow to relate the actual flow 
behind the detached shock wave to the subsonic (often 
incompressible) potential flow past the same body. 
pre- 


The 
dicted by several such theories is compared in Fig. | 
with the results of the numerical solutions described 
later. The standoff distance is a convenient test be- 
cause it is easily measured from schlieren photographs 


distance for spheres in air (y = 


and should be given accurately by any theory that is 
to predict other more important properties of the flow 
field. Moreover, its prediction is actually the main 
objective of most of these approximations. 


+ Another survey of the blunt-body problem in the special case 
of hypersonic flow will appear in the monograph Hypersonic Flow 
Theory by Wallace D 
published soon 


Haves and Ronald F. Probstein, to be 
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Fic. 1. Potential flow approximations for standoff distance of 


shock from sphere with 4+ 


The accuracy is seen to be mediocre, particularly in 
hypersonic flow (and is still worse for plane flow). 
Closest agreement is given by the theory of Kawamura, * 
who makes a straightforward match of streamline 
slopes behind the vertex of the shock wave with those 
His results have 
been reproduced by Heybey? who also attempts to esti- 
mate a correction for compressibility. However, it 
will be seen in section (1.3) that vorticity rather than 
compressibility is the chief factor neglected at high 
Mach Numbers and that Heybey’'s estimate of the com- 
pressibility effect is too large. In any these 
theories are of no value for predicting such quantities 


of the incompressible potential flow. 


case, 


as surface pressures (which by assumption have the 
form corresponding to incompressible potential flow). 


(1.2) Taylor Series Expansion From Shock 


If the form of a detached shock wave is known or 
assumed, the flow variables (pressure, density, velocity) 
just behind it are given by the oblique shock relations, 
and their first derivatives can be found by substituting 
into the equations of motion. Higher derivatives can 
be found by first differentiating the equations of motion. 
According to the Cauchy-Kowalewski theorem (see, e.g., 
reference 33) the flow field is analytic somewhat down- 
stream of the shock wave, and there is physically no 
reason to doubt that the region of analyticity extends 
to the body. One may therefore attempt to represent 
the flow layer between the shock and body by a few 
terms of a Taylor series expansion starting from the 
shock wave. This idea, introduced by Lin and Rubi- 
nov,’ has been applied to plane,* axisymmetric,®: 7!) 1! 
and three-dimensional’ flows. 

The expansion has been pursued furthest by 
Cabannes for axisymmetric flow. Taking the shock 
yave to be described in Cartesian coordinates by x = 
Ya,r", he gives the coefficients in the double Taylor 
series for the various flow variables out to 
(and of Stokes’ stream function to x‘r?, .. . ) for 
y = 7/5 and general free-stream Mach Number 1/ 
(reference 10), and the stream function to two addi- 
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tional terms in the special case .1/J = 2 (reference |] 


The variation of Stokes’ stream function along the ayjs 


of revolution behind a paraboloidal shock wave is com. ' 


pared in Fig. 2(a) with the accurate numerical solution 
Successive terms of the Taylor series are seen ty 
meander with no clear indication of either convergence 
or divergence. 

The expansion procedure can be simplified by intro. 
ducing the natural curvilinear coordinates associated 
with the shock wave 
oloidal shock wave, ete. 


parabolic coordinates for a parab. 
The Taylor series have then 
the advantage of being not double but only single ex. 
pansions (in the coordinate leading away from. the 
shock). 
gence properties of the series will show up more clearly 


It may also be anticipated that the conver- 


It is easy to recast Ca- 
The 


in the natural coordinates. 


bannes’ series in parabolic coordinates. result 


\\NUMERICAL 
10 .20 30 


(a) Cartesian coordinates. 


-.50 


.20 .30 


(b) Parabolic coordinates. 


Fic. 2. Variation of stream function along axis behind para- 
boloidal shock wave at JJ = 2 with 4 = 7/5. 
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corresponding to Fig. 2(a) is shown in Fig. 2(b), and 
it is now clear that at the nose of the body the series 
diverges. 

The probable source of the divergence has been sug- 
gested in reference 34. The Taylor series expansion 
describes not only the flow downstream of the shock 
wave but also its analytical continuation upstream. 
As indicated in Fig. 3, this fictitious flow will contain 
4 limiting line which intersects the axis at the sonic 
point and in the supersonic region is the envelope of 
outgoing characteristics. (Such a line exists also in the 
well-known Taylor-MeColl flow past a circular cone; 
its location is, in fact, tabulated by Kopal.*®) The 
flow variables are nonanalytic at the limiting line where 
they vary as half powers of the distance. Hence if the 
shock wave is closer to the limiting line than to the 
body, the Taylor series will not include the body in its 
radius of convergence. This suggestion has been veri- 
fied by running the numerical procedure backward to 
calculate the fictitious flow in the example of Fig. 2. 
On the axis, the limiting line was found to lie only 
three-fourths as far from the shock wave as does the 
body nose. 

The situation is comparably bad for other shapes and 
Mach Numbers and worse away from the axis or in 
plane flow. Hence the method of Taylor series expan- 
sion from the shock wave cannot be expected to suc- 
ceed. Nevertheless, the first three terms of the series 
(in natural coordinates) can serve as a useful qualitative 
model in the case of axisymmetric flow.** For a parab- 
oloidal shock wave at JJ = o (and y = 7 5) the 
model predicts the ratio of standoff distance to body 
nose radius correct to within 0.8 per cent, so that it 
may also have some quantitative value. However, it 
deteriorates away from the axis and does not predict a 
real body out to the sonic point. 


(1.3) Incompressible Approximation 


At high supersonic Mach Numbers the density varies 
only slightly between the shock wave and the body in 
the vicinity of the stagnation point. It therefore seems 
appropriate to solve the equations of rotational incom- 
pressible flow together with the exact conditions just 
behind the shock wave. Thus Lighthill'* has obtained 
the solution in closed form for the special case of a 
spherical shock wave at infinite Mach Number (that 
is, in the “‘strong shock’’ approximation), the body 
being found as a concentric sphere. Similarly, Whit- 
ham'* and Hayes” have shown that in plane flow a 
circular shock leads to a concentric circular body. 

To assess this approximation, numerical solutions 
have been carried out for a spherical and a circular 
shock wave in a stream of infinite Mach Number (all 
solutions for .1f = © were actually carried out at J = 
10,000) with y = 7 5. In axisymmetric flow the re- 
sulting b« nly shape, shown in Fig. 4(a), lies close to the 
sphere of the incompressible approximation out 
roughly halfway to the sonic point. In plane flow, as 
might be expected, the agreement is poorer (though it 


would improve in either case as y > 1). The body is 
roughly the front portion of a 4:1 ellipse standing on 
end rather than a circle. 

For the flow near the axis, these results were given 
earlier by Hida.'* His analysis is more general in that 
all supersonic Mach Numbers are considered. Fig. 
4(b) shows that his results for standoff distance are 
fairly accurate throughout the Mach Number range, 
the remaining discrepancy being presumably attribut- 
able to compressibility. Thus at infinite Mach Num- 
ber the approximation of irrotational incompressible 
flow*: * gives 0.094 as the standoff distance for a sphere 
in air; the increment due to rotation is 0.026 according 
to the Hida-Lighthill theory and the remaining differ- 
ence of 0.008 from the value 0.128 of the numerical 
solution must represent the effect of compressibility. 
(Heybey’s estimate of the compressibility effect’ is 
0.024 in this case.) It might be possible to include 
compressibility by iteration in the style oi the Janzen- 
Rayleigh approximation for subsonic flow. If so, this 
would seem the most promising of the analytic theories. 


(1.4) Newtonian Approximation and Improvements 
Thereon 


In simple Newtonian theory, fluid particles are im- 
agined to lose their normal component of momentum 
upon striking the surface of the body. The local pres- 
sure coefficient is then given by 2 sin* 6, where @ is the 
angle that the surface makes with the stream direction. 
Lees® has proposed a modified Newtonian theory which 
consists merely in scaling down this result so as to be 
exact at the stagnation point where the correct value is 
known. He suggests that because of compensating 
effects, this will be more accurate than the shock-layer 
or ‘‘Newtonian-plus-centrifugal” approximation'® 
which is the actual limit of the solution for .1/ ~ © and 


1. 
The surface pressures on a sphere and circular cylin- 
der at 1J = o~ with y = 75 according to these three 


appreximations are compared in Fig. 5 with numerical 
solutions. The modified Newtonian result is seen to 
fall either above or below the accurate solution, depend- 
ing on Mach Number and body shape, but to be re- 
markably accurate when y = 7 5. 

Chester?® and Freeman*! have taken the shock- 
layer or Newtonian-plus-centrifugal solution as a first 
step and sought to improve it systematically by suc- 
cessive approximations. The result is a double series 
expansion of the flow quantities for small 6 = (y — 1)+ 
(y + 1) and \/~? (a single series expansion in the work 
of Freeman, who considers only the ‘‘strong shock” 
case, corresponding to = ©). 

Unfortunately, the series appears to converge slowly 
for practical values of y, and worse in axisymmetric 
than plane flow. This is illustrated in Fig. 6 for the 
standoff distance in Chester's example of a parabo- 
loidal shock wave. By appeal to the model mentioned 
at the end of section (1.2), the present writer has sug- 
gested** that Chester's series actually converges for 
y > 2.2 (corresponding to 6 < 3 8), whereas for real 
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field upstream through shock wave, showing limiting line (from 
reference 34). 
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shock at M = ©, y = 7/5. (b) Standoff distance of shock from 
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SUPERSONIC BLU 


gases y does not exceed 1.67; but it is clearly not useful 
jor y above 1.1. In the same note it was shown how, 
by comparison with the model, the convergence of 
Chester's series for the standoff distance can appar- 
ently be accelerated. Thus at y = 1.667 in the ex- 
ample of Fig. 6, five terms of the modified series give 
AR, 0.198, which divers from the accurate numeri- 
cal result by only 4+ per cent. 
series is even. worse for surface pressure, and there no 


The convergence of the 
useful modification has been found. 


(2) DESCRIPTION OF NUMERICAL METHOD 


The preceding survey shows that none of the existing 
analytic treatments is adequate for calculating the de- 
One must turn instead 

A rela- 
described 


tails of flow near a blunt nose. 
to numerical solution of the full equations. 
tively simple numerical scheme will be 
which has been programed for a medium-sized elec- 


tronic computer. 


(2.1) Coordinate System 


It is convenient to assume a family of detached shock 
waves of known shape. The resulting simplification 
more than offsets the disadvantage of having to accept 
whatever body shapes result. Another possible objec- 
tion will be dealt with in section (3.4). Schlieren pictures 
show that the shocks produced by simple shapes such 
as spheres, paraboloids, and ellipsoids (or their plane 
counterparts) are themselves nearly conic sections. 
Accordingly, the present method has been applied to 
the family of plane and axisymmetric bodies that 
support detached shock waves described by conic 
hyperbola, parabola, prolate ellipse, circle, 
It could, however, be extended to 


sections 
or oblate ellipse. 
other analytic shock shapes. 
In Cartesian coordinates originating from its vertex 
(Fig. 7), any conie section is described by 
r? = 2Rx — Bx? (1) 
Here KX is the nose radius, and B the bluntness, which 
is a convenient parameter that characterizes the eccen- 
(The bluntness B is a? 
As shown 


tricity of the conic section. 
for an ellipse but —b? a? for a hyperbola.) 
in Fig. 8, the bluntness is zero for a parabola, negative 
for hyperbolas, positive for ellipses, and unity for a 
circle. 

It is advantageous to choose a natural coordinate 
system that contains the shock wave as one of its sur- 
faces. The conventional hyperbolic, parabolic, elliptic, 
and polar coordinates are unsatisfactory because their 
definition changes at B = 0 and again at B = | whereas 
conic sections clearly form a single continuously vary- 
ing family. A suitable unified orthogonal system that 
covers the entire range is introduced by setting 


r = 


Il 


l'he shock wave of bluntness B, is described by 7 = 1; 
the downstream axis by 7 = 0; and the upstream axis 


°T-BODY PROBLEM 
by = 0 (Fig. 7). Special cases are 
vr R= (1 2)(1+ & — n°), for parabola (B, = 0) 
V1 — for circle (B, = 1) 


For ellipses, Eq. (2) gives only the left Fall which is the 
part that may form a shock wave. The right half of 
the ellipse is obtained by taking a plus instead of a 
minus sign ahead of the radical in Eq. (2). 

In this coordinate system the line element 1s given by 


= hydt? + he?dyn? + (3a) 
where vy = 0 for plane flow and vy = | for axisymmetric 
flow, so that the last term (involving the azimuthal 
angle ¢) appears only in the latter case, and 
hy? = (Ce + 9?)/(1 — BE), 

he? = (CE + ?)/(C + Bon’), hs? 


— 


(3b) 
where C= 1 


(2.2) Equations of Motion 

Let |’ be the velocity referred to the free-stream 
speed |”... p the density referred to its free-stream value 
p., and p the pressure referred to p.. 1 Then the 


r r2= 2Rx-Bx2 
€=0 
a 
Fic. 7. Coordinates for conic section. 
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Fic. 8. Conic sections of varying bluntness. 
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[(En)” V (C2 + (C + Bn*)pu', + 


[(&n)’ (CE + n°) (1 — Bé*)pv], = 0 
STAGNATION POINT p( ur; — [Cév? (Ce + + 
WITHOUT V(C + Bn’), (1 — Be) + 
ITERATION 
6.3 / WITH EXTERNAL Inu + + = 0 
ITERATION (5 
— [nu? (Ce + n°)] + 
Vil Bé&) (C + Bn?) + 
[Cév (Ce + +p, = 0 
100 96 .92 88 .84 80 76 72 
7 vV C + Bn? (p p’), = 0 ) (] 
Fic. 9. Variation of density along axis behind parabolic shock = = 
wave at M = ~, y = 7/5. where uw, v are the components of |’ in the &, 7 directions Cs 
The first (continuity) equation is satisfied by intro- | 
equations of continuity, motion, and energy are, in ducing the stream function V according to wh 
vector form, y, 
Vv, = (én)’ V (CE + 9°) (C+ Bn*)pu | 
div(pl’) = 0, p(I’-grad) + grad p = 0, (6 
¥, = — (Ce* + 97) /(1 — 
V-grad(p = 0 (4) 
Phen the last (energy) equation simply states that : bee 
where y is the adiabatic exponent. Transforming the 
= (,) 
these to the (£, coordinate system (using Eq. (3), and Ovi 
dropping the subscript on B, for simplicity) gives Using this to eliminate the pressure from the equa- | _ line 
tions of motion gives } eva 
pro 
— — p) + — Vi + (v — (CE + (WE + 
[((C + Bn®)/ (1 — + [BEC — Be) Wav, + 
[(CE? + — Bé*)] fo, + = 0 | ( 
— Bé) (C + Bn?) + — [Bn (C + Bn?) + val 
[(CE® + + + = 0 
( 
(2.3) Initial Conditions 
Values of a, v, p, p just behind the shock wave (at po = (y + 1) — Be) [21 + Ce) + 
n = 1) are found from the oblique shock relations | 
e.g., reference 37—in terms of the slope of the bow (y — DAP — Bé)], (9) tint 
wave, V1 Be Expressed in terms of the stream = + »), v, = aty = 1 
function, these give the initial conditions p {a 
put 
and for the function tery 
by 
2 
2y 1/71 — Bs?) — (y — + Cs?) 
(v) = x p 


y(y + (1 + Cs?) 


[2(1 + Cs*) + (y — — | 
(y + — Bs?) [1 + »)¥] | Fig, 


fron 

(2.4) Form of Problem for Numerical Computation T 
i it i tati 

For numerical work it is advantageous to use (1 + v)W ¢'*", which is constant on the shock and elsewhere more a 
nearly independent of ¢ than is VY and vanishes only on the body rather than also on the axis of symmetry. Hence set the 

n) = + 0) (11) spac 

Stra 


Then the initial conditions become the 
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and the equations of motion 
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aty = 1 


+ 
3 (w+ 


at 
(13a) 
# 


l+p n C + Bn? l+p pp 


> 


C2 + 7° C+ l+yp l+yp 


where f= 


y(y + 1)M? (1 + Cs?) 


(2.5) Numerical Integration 


The initial value problem of Eqs. (12) and (13) has 
been solved numerically by forward integration from 
the shock wave (y = 1) toward smaller values of 7. 
Over each interval (An), p and w, are extrapolated 
At each new value of n the £-derivatives are 
The 


linearly. 
evaluated by 11-point numerical differentiation.® 
procedure may be summarized as follows: 

(0) Calculate initial values at » = 1 from Eq. (12). 

(1) Calculate &-derivatives pz, w:, wy, by 11-point 
numerical differentiation. 

(2a) Calculate p, from Eq. (13a). 

(2b) Calculate w,, from Eq. (13b). 

(3a) Extrapolate p and w, linearly to next smaller 
value of 7: 
m+1) 


p” +1) — p (An) = (An), 


(Sb) Extrapolate w using averaged value of w,: 
= — (1/2) (ag) + 

(4) Repeat steps (1) to (4) at new value of » and con- 
tinue until w is negative. 

After completion of the integration, the accuracy of 
p (and hence of p) can be greatly increased by recom- 
puting it using the averaged value of p, over each in- 
terval \y, so that the value after the mth step is given 
by 


— (An) 2)p,“ + a” + + 


= p® 
(14) 


Fig. 9 shows in a typical case the improvement obtained 
from this ‘‘external iteration.”’ 

This procedure has been coded for machine compu- 
tation on the IBM 650 electronic digital computer with 
floating decimal attachment. For each value of n 
the flow quantities are calculated at the 20 equally 
spaced values £ = (nw + 1 2) (AE), O < wm < 19, which 
straddle the axis of symmetry £ = 0 in order to avoid 


the indeterminate form p: &p in Eq. (13). The mesh 


+ (1 + fy? 


+ 


(13c) 
(y + 1).1/°(1 — Bs?) 


dimensions Aé and An are arbitrary (but equal for all 


steps). One prescribes the parameters 

M = free-stream Mach Number 

B, = bluntness of conic section describing shock 

v = (for plane, | for axisymmetric flow 

7 = adiabatic exponent 

Aé, An = mesh dimensions 
Flow variables (and their n-derivatives) are printed for 
each mesh point. The body shape and surface pres- 
sures can then be found by interpolation. Machine 
computing time is 1! , minutes per value of yn. It will 


be seen that + to 8 steps in 7 yield ample accuracy so 
that a typical case requires 5 to 10 minutes computing 


time. 


(2.6) Instability 


In the subsonic region the initial value problem is 
unstable. This means that small errors (such as arise 
from rounding to the 8 significant figures used) grow 
in geometric progression. The instability has been 
tested in a typical case by perturbing the initial data, 
changing w from 1.000 0000 to 1.000 0010 at the point on 
the shock wave nearest to the axis (where the equation 
is most elliptic and hence most unstable). Fig. 10 
shows that in 11 steps (last step interpolated) from 
shock to body a unit error grows by a factor of 160. 
This means that rounding errors of | 2 in the last 
place can amplify to invalidate no more than 2 of the 
S significant figures used. Near the stagnation point 
the error grows by a factor of about 3 in each step. 
(This may be compared with the maximum factor of 
5.8284 for Laplace's equation in a_ semi-infinite 
domain.**) Hence the error would swamp all 8 signifi- 
cant figures if it were necessary to take as many as 24 
steps from shock to body. Fortunately, many fewer 
steps than even 11 suffice, so that instability is of no 


practical concern. 


fw Bs 7 
- 
1+ — p 1+ yp 
|_| fw, \? | C+ Br*/ ke, \* 
+9 l+~y» 1 — Be\l + 
: 
1 — Be 
(7 
= | aa 
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Another type of instability peculiar to the particular 
differentiation scheme used is more serious and, in 
practice, limits the downstream extent of flow field that 
can be calculated. The 1|1-point differentiation scheme 
uses central differences where possible, but the largest 
5 values of € must use progressively more noncentral 
schemes. These involve large differences of small num- 
bers and so produce instability which is more severe in 
the hyperbolic than the elliptic region. It eventually 
causes oscillations of flow variables for the highest few 
values of &, and stops the computation by producing 
negative values of p. This difficulty has been allevi- 
ated by using 7- and 5-point differentiation schemes at 
the highest two values of &. It is to be entirely elimi- 
nated by recoding the program to use only central-dif- 
ference schemes, starting with more values of € and 
dropping the highest 5 in each step in 7. 

A second imperfection in the present procedure ap- 
pears when the factor }w + (fw; (1 + v)]} in Eq. (13) 
vanishes, leading to a singularity in w,,. This can hap- 
pen only inside the body (for w < 0) but, because of the 
instability illustrated in Fig. 10, its effect occasionally 
spreads so rapidly as to affect the actual flow field. 


(2.7) Accuracy Comparison With Garabedian- 

Lieberstein Method 

The accuracy of the numerical method has been 
tested both internally, by refining the mesh size, and 
externally, by comparing with experiment and with the 
numerical procedure of Garabedian and Lieberstein.” 

The Garabedian-Lieberstein method is like that pro- 
posed here in starting from a given analytic shock shape 
(and conic sections were actually used) and solving the 
initial-value problem with electronic computers. It 
differs in that it avoids the instability in the subsonic 
region by continuing the initial data analytically into a 
fictitious third dimension where the equation is hyper- 
bolic rather than elliptic. The solution is carried out 
by the numerical method of characteristics on a num- 
ber of planes in the fictitious three-dimensional space 
that intersect the real physical plane along curved 
lines. Thanks to known uniqueness proofs for the nu- 
merical method of characteristics, the Garabedian- 
Lieberstein method possesses a degree of mathematical 
rigor seldom attained in aerodynamic research. 

The price paid for eliminating instability is a large 
increase in the computing time (and cost) compared 
with the present method. Also the solution is arti- 
ficially restricted to the subsonic region (though an- 
other solution can be carried out for the supersonic 
region). 

It can be shown that the present procedure converges 
as the mesh size is refined, despite its instability. In- 
deed, it differs only in degree from any stable numerical 
procedure, in which one must strike a balance between 
the truncation error resulting from too few steps and 
the round-off error resulting from too many; here the 
accumulated error simply grows in geometric rather 
than arithmetic progression. What is lacking is a proof 
that the solution to which the present method con- 
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verges is the correct one. (The question of whether 

a real body exists for any given shock wave is an oper 

question in either method.) However, it reproduces ! 
the results of the Garabedian-Lieberstein method t 

four significant figures, as shown for one case in Figs 

11 and 12. Fig. 12 indicates that the present method 

is the more accurate near the sonic line. 

In general, as in the example of Fig. 11, four to siy 
steps in 7 along the axis are found to yield the standof 
distance correct to within | per cent, with comparable 
accuracy in the other flow quantities. It was thought 
worthwhile to carry out the solutions discussed beloy 
only to this standard of accuracy. 


(3) EXAMPLES AND DISCUSSION 


Some 50 examples have been calculated, covering « 
wide range of shock bluntness, Mach Number, and ; | 
in both plane and axisymmetric flow. Only axisym. 
metric flow will be considered henceforth as being oj 
more practical interest. A complete compilation of both 
plane and axisymmetric solutions will be issued later 
as an NACA publication by the writer and Helen 
Gordon, who carried out the coding and machine com- 


putation.* 


(3.1) The Family of Bodies 


At each Mach Number (and y) the family of conic- 
section shock waves leads to a one-parameter family of | 
bodies. Fig. 13 shows various members of the family 
for JJ = ©. The downstream extent of the body was 
limited in each case by the end instability discussed in 
section (2.6), which could be eliminated. 

Except possibly for the very bluntest, the bodies can | 
all be accurately described by conic sections back at 
least to the sonic point or limiting characteristic. 
Thereafter the solution could be continued, or modi- 
fied, by the method of characteristics. However, it 1s 
a virtue of the present method that it gives the sub- 
sequent shape. The sonic line and the change of type 
from elliptic to hyperbolic play no significant role 
From this point of view it may be said that the blunt- 
body problem is not a transonic problem. 


(3.2) Comparison With Experiments on Spheres 
Because the bodies are found all to be represented 

closely by conic sections, it is possible to consider é 

fixed shape over a range of Mach Numbers. Enough 


experience has been accumulated at y = 7/5 that a de- 
sired body can be produced in one or two attempts. 
Thus the sphere has been solved for a number of Mach 
Numbers from 1.30 to © (and values closer to 1 would 
offer even less difficulty). The standoff distance 1s 
compared in Fig. 14 with experimental results from 
references 40 to 44. The surface pressure distribution 
is compared in Fig. 15 with experiments at / = 1.5 
(reference 45) and ./ = 5.8 (reference 41). 


* We are indebted to Marcelline Chartz for continued guidance 
through the mysteries of machine computation. 
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Variation of body bluntness with shock wave blunt- 
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Fic. 16. 


(3.3) Comparison With Other Numerical Methods 


After developing the present numerical procedure, 
the author learned of two other independent treatments 
by essentially the same method. Priority clearly goes 
to K. Mangler, of the Royal Aircraft Establishment, 
although his work is not yet published.* At a meeting 
in Holland in June 1957, he reported a hand computa- 
tion of the plane flow behind a parabolic shock wave 
at lJ = 7 (with y = 7 5). The body was found to be 
circular within half a per cent to behind the sonic line 
(in conformity with results of the present method) 
and was continued to a semicircle by the numerical 
method of characteristics. 

More recently, Zlotnick and Newman* have applied 
a similar procedure to the spherical shock wave at 
several high Mach Numbers. The accuracy of their 
examples is low (the pressure at the stagnation point 
being off by some 3 per cent and the standoff distance 
low by 10 per cent), at least partly because the thick- 
ness of the layer between shock and body is neglected 
at one point. They apparently attempt to suppress 
the instability by filtering out high harmonics (which 
are the most unstable). The present author believes 
such attempts to be dangerous as well as unnecessary. 

In contrast with these indirect methods (that start 
from an assumed shock wave), Belotserkovsky*! and 
Uchida and Yasuhara**? have described methods for 
solving the direct problem of a given body. Following 
the so-called integral method of Dorodnitsyn, Belot- 
serkovsky divides the layer between shock and body 
into # strips (2 or 3 in practice) across each of which 
the flow variables are approximated by polynomials. 
Integrating the equations of motion analytically across 


* ADDENDUM: Mangler’s work is now available in: Mangler, 
kK. W., and Evans, M. E., The Calculation of the Inviscid Flow 
Between a Detached Bow Wave and a Body, R.A.E. Tech. Note 
Aero 2536, October, 1957. 
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the various strips yields a system of ordinary differ. | 
ential equations that are integrated on an electronj 
digital computer. As in the Garabedian-Liebersteiy | 
method, the sonic line assumes a spurious significance 
which necessitates special treatment in its vicinity 
Belotserkovsky has published results for a circular 
evlinder at JJ = 3, 4, 5 (with y = 7 5) which agree 
well with those from the present method; and it would 
seem reasonable to assume that he has also considered 
more practical axisymmetric bodies. 

Uchida and Yasuhara have given a very tedious 
method of successive approximations which does not 
seem suited to machine computation. Their one pub. 
lished example of a circular cylinder at J = 2(y =7 5 
agrees fairly well with results from the present method, 


— 


(3.4) Shock Insensitivity 


An objection sometimes raised against the inverse 
method of starting from an assumed shock wave is that 
its shape is insensitive to changes in the body shape. 
A near nonuniqueness is implied according to which 
essentially the same shock wave would correspond to 
However, no such difficulty exists J 
Fig. 16 


— 


various bodies. 
within the family of shapes considered here. 
shows that for shapes sharper than a parabola the | 
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with same nose radius and with same sonic point : 

Fic. 18. Result of trying constant effective y to simulate real 

gas effects, = 14.2. 
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bluntness of shock and body change at the same rate; 
and although the shock insensitivity increases for 
blunter shapes, it by no means approaches infinite slope 
which would imply complete insensitivity. 


(3.5) Shoulder Choking 

Increasing shock insensitivity for blunt bodies is 
related to the phenomenon of shoulder choking pointed 
out by Busemann** and Hayes.'* According to this, 
the standoff distance and the shock shape in the sub- 
sonic region are determined mainly by the shape 
of the body near the sonic point. This is illustrated in 
Fig. 17 by superposing first the vertices and then the 
sonic points of two different bodies. The effect would 
become even more pronounced for blunter shapes (and 


in plane flow). 


(3.6) Possibility of Treating Real Gases 


The present numerical method is suitable for the in- 
clusion of real gas effects assuming equilibrium thermo- 
dynamics. Hand calculation would probably be neces- 
sary if the standard tables were to be used, but they 
could be approximated analytically for purposes of 
machine computation. 

Several investigators have suggested the use of a con- 
stant effective y to account roughly for real gas effects. 
This idea was tested by starting from the observed 
shock for a hemisphere flying at .1/ = 14.2 (reference 
7), which clearly involves significant gas imperfec- 
tions, and attempting to reproduce the body. The 
observed standoff distance is obtained by choosing y = 
1.35 which according to section (3.5) must be the ef- 
fective value in the vicinity of the sonic line. Fig. 1S 
shows that the body shape is then closely reproduced. 
At least in this case of mild gas imperfections, the idea 
shows promise. 
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The Flow of a Supersonic Jet in a Supersonic 
Stream at an Angle of Attack 


F. EDWARD EHLERS* anp TORSTEIN STRAND ** 


Boeing Airplane Company 


SUMMARY 


The lisearized equations for the flow induced by a supersonic 
jet at a small angle of attack in a supersonic stream have been 
solved using the Laplace-transform theory. Reflections from 


the jet interface depend upon a parameter 
k= V (Me? — 1)/(M? 1) 


where ./ is the Mach Number, ? the pressure, and the sub- 
scripts, | and 2, refer to the jet stream and the exterior flow, 
respectively If k = 1, for two dimensions there are no reflec- 
tions from the jet interface and the jet aligns itself with the ori- 
ginal undisturbed flow in a finite distance. For k > 1, the jet 
interface for the symmetrical solution oscillates, but the slope 
is monotonically decreasing for k < 1 However, the angle of 
attack term shows an oscillating boundary only for k < 1. 

The three-dimensional solution has no oscillations of the 
boundary ifk = 1. For other values of &, the sum of the residues 
of the Laplace transform introduces singularities in the solution 
which appear to be a result of the linearization. Simple relations 
are obtained for the perturbation velocity potentials near the 
Mach cone from the rim, x — B.(r — 1) = O, and far downstream 
of the jet exhaust. The pressure discontinuity decays like 
1/Vr along the Mach line from the rim of the jet. 

The length of the jet engine in three-dimensional flow influences 
the behavior of the jet. For a tube length of infinity, the jet 
far downstream behaves like an infinite cylinder at an angle of 
attack; but for a tube length of zero, the angle of the jet far 
downstream is equal to a/(1 + 6,), where a is the angle of attack 
of the undisturbed exterior flow and 6; is the ratio of dynamic 
pressures of the jet stream and the interior flow, respectively. 


SYMBOLS 


P = static pressure (with subscript = initial static 
pressure ) 

M = initial Mach Number 

8 = 

V = velocity anywhere in flow where pressure = P? 

U = initial velocity 

a = angle of attack of undisturbed flow with respect 
to jet-tube axis 

Ap = angle of flow deflection at r = 1 

x,”,8 = cylindrical coordinates (@ = zero along top of jet) 

x,y = Cartesian coordinates 

to) = total velocity potential 

= velocity-perturbation potentials 

6 = ratio of dynamic pressures = ¥;.7,?P;/y2Me?P» 

€ = (P, — P2)/y2P2M? 

p = density 

k = 

p, A, s = Laplace-transform variables 


Received June 6, 1955. Revised and received December 7, 
1956. 
* Research Specialist, Boeing Scientific Research Laboratories. 
** Formerly, Aerodynamicist. Now, Senior Aerodynamics 


Engineer, Convair (San Diego), A Division of General Dynamics 
Corp. 


| e » (x)dx 
J0 


I, (pB) = modified Bessel function of the first kind of order v 
K,(p8) = modified Bessel function of the second kind of 


order v 

Dip) = + )Ki( pps) 

Exp) = asymptotic form of D.(p) 

h = small Mach Number parameter = 6; — By 

¥¢o = Ward's solution for flow over tube at angle of attack 

go” = additional velocity-perturbation potential due to 
jet for exterior flow at angle of attack = go’ — gp 

K,E = elliptic integrals of the first and second kind, re- 
spectively 

= x/28 

ne? = |x — B(r — 1)]/28 

= (1 —k)/(1 +h) 


eNK if Ap, 
= 
J Br, MA — 


1 
dy 
Br, K,"(AB2) 


Subscripts 1 and 2 denote quantities of the jet interior and 
exterior flow, respectively (exceptions 6, f, g, and the Bessel func- 


V(1/B2) = 


tions). 
Superscripts (0) and (1) denote the symmetric and nonsym- 
metric solutions, respectively. The parentheses and the 0 are 


dropped when used with ¢. 


INTRODUCTION 


| Bs US REPRESENT a supersonic jet engine in super- 
sonic external flow by a hollow circular cylinder of 
radius ry = | and length / at an angle of attack a (see 
Fig. 1). Furthermore, this angle of attack and the 
difference in pressure between the two streams, (?,; — 
P:) are assumed small, and hence the linearized theory 
applies. The suction at the jet-engine intake must be 
such that there is no spill-over at the entrance. The 
flows inside and outside the tube then may be treated 
separately. Issuing from the exit of the jet is a uni- 
form and parallel flow of Mach Number .1/; and pres- 
sure /;, while in front of the cylinder there is a uniform 
parallel flow of Mach Number .1/, and pressure 7. 

The method of solution is to simplify the linearized 
boundary conditions and the partial differential equa- 
tions for the perturbation potentials by means of the 
Laplace transformation. The Laplace transforms of 
the potentials then become solutions of ordinary second- 
order differential equations. 

For completeness, the two-dimensional solution for 
the inclined jet is also given. The results for zero angle 
of attack are similar to that of Pai? and of Pack.* 
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Fic. 1. Flow model and coordinate system (three-dimensional). 


THEORY 


Let us assume that the velocity potentials for both 
the jet and the exterior stream can be expressed in an 
expansion of the form 


n=0 
where LU is the uniform stream velocity, a the angle of 
attack of the jet tube, x, 7, @ are the usual cylindrical 
coordinates, and y‘"’ the perturbation potentials. We 
shall use subscripts | and 2 to denote the quantities for 
the jet interior flow and for the exterior flow, respec- 
tively. To satisfy the continuity of pressure at the 
jet interface, we shall require the expansion of (1° + 


U*) — 1 in powers of a. Thus, 


n=0 
@o 


= + Day}? + 


n=0 n=0 


Expanding and simplifying leads to 


(V2 U?) + + 
where 
QO = 2y, 02 (Wo? y,? 
2y,") + 2y,y, + 
+ (2 WoO Po? 
= + WY +p? + +... 


Now the pressure for the jet flow is given by 


which near the boundary may be expanded in the form 


P= Pi + (V2 + 
8) Pill — +... 


Substituting for 1 — (1° U4") and writing the corre- 
sponding quantity for the exterior flow, we find that 
the condition requiring that the pressures be equal on 
the jet interface becomes 
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P; ("1 2)M [Qi + 4+ + 
Py — (y2M2?P2/2) + aQe +...) 4 
+ 2aQe Qe +... 


This can easily be written in the form 


— — 

Q2 + aQ +... 4+ 

(Qi? + +...) — (Me?/4) 

+ 4+... 


X 
X 


Letting 


= 6, (Pi — y2P2M? = «¢ 
= My? 4 = 


Since this must be true for all a, then we must set the 
coefficient of each power of a equal to zero. This 
yields a sequence of equations, of which the first two 
are 


+ + 26.0,°Q — = 0 


If we neglect all products and powers of the potentials 
y and their derivatives higher than unity, then the two 
equations above reduce to 


Now the condition that the jet and the outside stream 
have the same inclination at the interface is expressed 
by 


(1) 


= P2, P», 


on the common boundary. This becomes 


+ apis’? por + ayo +... 
Multiplying the equation by a product of the denomi- 


nators and neglecting products of the perturbation po- 
tentials yields 


Wi,” + ay,,\) = Yo, 


Equating coefficients of the powers of a@ leads to the 
following boundary conditions for the first two approxi- 
mations: 


Eqs. (1) and (2) are the boundary conditions for the 

jet interface to the first order in a. Since we have 

assumed the disturbance of the jet to be small, we shall 

satisfy the boundary conditions at r = 1, the mean 

position of the jet boundary. 

For the three-dimensional jet, we let 

= cos O[r + 
(3) 


= cos 0¢;'(x, 7), 


The term r cos 6 of Y2"" is the potential for the upstream 
components of velocity due to the angle of attack, and 
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g./(x, 7) is the perturbation caused by the inclination 

of the jet. With these substitutions, the boundary con- 
ditions, Eqs. (1) and (3), become 

€ + G2 = O1¢1r (4a) 

Gor" (4b) 

and = $2 (5a) 

= l + (Sb) 


The development for the two-dimensional jet follows 
similarly with 


W(x, 7,8) = (x, r, 0) = 


W(x, 7,8) = 
y+ ¥) 


The boundary equations are the same as Eqs. (4) and 
(5) with y replacing r. 

Since the perturbations of the two streams are as- 
sumed to be small, the boundary conditions for the 
three-dimensional jet are satisfied at r = 1, the mean 
location of the jet interface, and at y = +1 for the 
two-dimensional flow. 


THE AXIALLY SYMMETRIC THREE-DIMENSIONAL FLOW 


The axially symmetric solution has been found by 
D. C. Pack in reference 6, and the reader is referred to 
that paper for a complete discussion of results — In our 
notation, the potentials are given by 

Boe pBir) Ki - 

‘ 
J ,-; p?D\(p) 

€ 
p?D\(p) 


where 


Di(p) = 11(pB1)Ko( + Ki( (9) 
k = (62/8:)6; and B = VM? — 1 


Simple expressions for g; and ¢g2 in the neighborhood of 
the Mach cones x + 8:(7 — 1) = Oand x — 6.(r — 1) = 
0 can be obtained by using the asymptotic forms of the 
transforms. For large s we have 


Ki(t)~ Va 22e, 


(e° + Vs) 
Be 
pl + 


Be lx + — 1)] 


Similarly, 
elx — — 1)] 


A SUPERSONIC JET 499 


At the rim of the jet, x = 0, 7 = 1, we have 


Or|,=1 = Boe (1 + k) (10) 


FIRST-ORDER APPROXIMATION 


As a first-order approximation,* we assume that the 
differences in both the pressure and the Mach Number 
are small—i.e., 


(Pi — = Bi — =h 


where e, 4 are small quantities. In this case, 
k ] + Ole, h), = Ve 


It should also be pointed out that this expression is 
the exact solution for 


k=1 Bi = Be 


Then the flow issuing from the tube will have the same 
Mach Number, and the pressure and ratio of specific 
heats are not equal but must satisfy 


and 


= 


Let s = pp, then, expanding the Laplace transform of 
¢2 in powers of « and / and neglecting higher-order 
terms than the first degree leads to 


f (11) 


Substituting 


T(s) = 


l scos@ 
€ cos 6d6 
TJ0 


interchanging the order of integration, and evaluating 
the transform of Ao(s7), we get 


= X 


or cos’ '(—x/ B+r) 
f [sin® @ V (x B+ cos 6)? — 
0 
(12) 


with upper limit 


x — B(r 1) > 28 
x — Bir — 1) < 28 


= for 
= cos''(—x r) for 


Similarly, 


= (8 m) X 


x or cos”! [—(x/B+r)/r] 
[V(x 8+ rcos 6)? — 1] (13) 
0 


It is noted that, to this order of approximation, the 
pressure difference determines the flow, and the differ- 
ence in Mach Numbers enters only as a second-order 
correction. 

In order to find the shape of the jet, Eq. (11) was 
differentiated with respect tor. It can be shown that 
for r = 1, ¢, may be represented by elliptic functions 
in the notation of reference 7. The velocity component 

* The first-order approximation employed here was suggested 
by Prof. J. D. Cole, California Institute of Technology. 
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Fic. 2. Contours of integration. 


¢g», from Eq. (11) on the boundary of the jet becomes 


r=1 


Now from Watson (reference 4, p. 44, Eq. 4), we have 


T(s)Ky(s) = —(2/m) Ko(2s cos @)cos 26 dé 
0 


Substituting into g., and interchanging the order of 
integration, we get 


= X 


|cos 26/V (x 28)? — cos? 6] for x/28> 1 


r=1 
{ |cos 26 V (x 28)? — cos? 0]d6, for x/2B < 1 
cos ~'(x/ 
To convert the first equation to elliptic integrals for 
x ‘28 > 1, we write 


= (68m) X 
[(2 sin? — 1) mV 1 — (1/m)? sin? 
0 


where m7, = x/28. By elementary algebraic manipu- 


lation this becomes 
gor = (€B/m) [(2m? — 1)K(1/m)/m — 2mE(1/m)]) 


For x/28 < 1, we let cos 6 = n; sin ¢; then 


= (€B/m) x 


n/2 
f [(2m? sin? — 1)/W1 — m? sin? 
0 
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or go, = [K(m) — 2E(m)] 


where 9; = x/2¢. 

The above relations have been plotted in Fig. 3. By 
integration, the streamline shape and location were ob. 
tained and are shown in Fig. +. At x« = 2¢ the deriy- 
ative g», becomes infinite. In fact, this singularity js 
seen from Eqs. (12) and (13) to occur on the character. 
istic x/8 — (ry — 1) = 2 for both the exterior and in- 
terior flows. This is the Mach cone emanating from 
the point on the axis x = 6 where the disturbance from 
the rim of the jet is focused. These Mach lines appear 
to be most predominant in photographs of a supersonic 
axially symmetric jet in a supersonic stream (see refer- 
ence 6). 


THE ANGLE OF ATTACK SOLUTION FOR THE THREE- 
DIMENSIONAL JET 

When the potentials of Eq. (3) are substituted into 
the differential equation for y, the differential equations 
for and take the form 

The solution of this differential equation for the flow 
over the cylinder has already been given by Ward.' 
For our particular situation (see Fig. 1) it takes the 
form 


A(x4+l) 
e Ky (Aer) 
fry (AB») 
This satisfies the boundary condition Og) Or = —|I at 


r=lfor—/<vx. If welet 


a | | 
| | 
aja | 
rolw 
| \ 
| 
| | | 
| | 
° 
—- 
P 
-.2 
Fic. 3. Slope of streamline = yg, atr = 1, a = 0. First-order 
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then the boundary conditions, Eqs. (4b) and (5b), for 
the Laplace transforms become, forr = 1, 


= So — [¢o(0)/p] + 


and Gir = $2,” 
Since, atr = 1, 
0 Bri (AB2) f 
a 
J Br Jo { 
(AB2)dd 
then 
ao — p] = p (14) 
where 


V(p) = 


2Qri 


Now the differential equations for the Laplace trans- 
forms of ¢g;’ and ¢»” satisfy, respectively, 


Sore” + (L/r) ger” — + (1/p?)]e2” = 0 

The general solutions for 2)’ and g»” then take the forms 
= folp) (phir), = go(p)Ki( 


since ¢;’ must be finite at y = 0 and ¢g.” must represent 
outgoing waves at infinity. Substitution into the 
boundary conditions yields 


(82 81) (pB2)go(p) — = 
Ki(pB2)g(p) — = Pp] 


These equations enable us to determine /, and go, 
which, after substitution and inversion, yield 


A (15) 
J Br: pD2(p) 
where 


Since V(p) is the transform of O¢go/Ox'|,=}, it is possible 
by using the convolution theorem to write g»” in the 
following form 


0 


ex —Be(r—1) 
| F(x — 1) V(t/Be)dt 


where 


vie) = eK 


ae + 
> 
im 
< 
& - 
q 
¢ e ‘2 
« - 
2p 


Fic. 4. Streamline atr = 1,@=0. First-order approximation. 


is Ward's solution for the pressure on the outside of the 
tube and 


F(x —t)= dp 


| p81) Ki (pBer) 
Br, 

When the Br, contour is replaced by Br. contour (see 
Fig. 2), periodic terms result from the residues at the 
poles of the integrand given by the complex roots of 
the equation D.(p) = 0. 

The function F takes the form 

r{l + (61 B»)R] 

0 \[A?+ 2° 


F(z) = §)(z) + 


where 
Ty — 


A 
B 


and S,(s) represents the sum of the residues at the poles 
of the transform given by the complex roots of D2(p) = 
0. For large p, the Bessel functions may be approxi- 
mated by the first terms of their asymptotic expansions. 


Thus, 


aa” 
ef = 


(27 pei 2 


Il 


= Ex(p) 


where the upper signs hold for the upper half plane and 
the lower signs for the lower half plane. The roots of 
this approximate relation are 


= (1/2) log + — (1/4) sgn (1 — 


where wu = (1 — k) (1 + k). These roots are similar 
to those given by Pack® for the symmetrical solution. 
However, we note that the imaginary part is 


— (1/4) sgn (1 — 
for the angle of attack potential and 


x + (1 4) sgn (1 — 
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for the symmetrical solution. By an argument almost 
identical to that given by Pack on pages 750-751 of 
reference 6, it can be shown that there is a one to one 
correspondence between the roots of D.(p) and of its 
asymptotic form /,(p). It becomes immediately ap- 
parent that there are no poles in a finite region of the p 
plane if k = 1. 

Replacing the Bessel functions by their asymptotic 
forms, we get, for the Laplace transform in the upper 
half p plane, 


(1 + k) Vrp(e?® + in) 
The poles are at 
or = (1/2) log + — (1/4)] 
if « > 0. The sum of the residues, approximated by 
their asymptotic forms, becomes 
—2k 
(1 — 
This may be written 
2k 


— B2(r — 1) 
f Real ( dt 
bi 


where the contour is chosen for real z — B.(r — 1) > O 
so that the integral and the sum converge. Substi- 
tuting for p, yields 
2k 

x 

V — 

f el! 28)) log Real (int/4p,) dt 


n=1 


Real > 


So = 
n=1 Pn 


So(z) = 


So(z) & 


Summing the series and simplifying yields 
k 


x 
— k?) 
sin (mt/2) 


So(z) 
dt 


The function has logarithmic singularities when 
sin (7/2) = O and sin (32t/4) 0 
The points of singularity are 


[z — B(r — 1)]/B = 4m — 2. whenm = 


Similarly, the velocity potential for the interior flow 
becomes 


x+Bi(r—1) 
| G(x — t)V(t/B.)dt 
0 


where 
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r 
[1 + (81 


J, 
0 + 


and S;(z) is the sum of the residues of the transform 


Bol (pBir) Ky’ (pBe) 


17 
Bi pD2(p) 


at the roots of D.(p) = 0. The nature of the singulari- 
ties of the interior solution are found in the same way 
as for the exterior solution. Replacing the asymptotic 
forms of /; and A, for the upper half plane into Eq. (17), 
we get 


Bi pr/r (1 + + ig) 
The residues are 


j 


Real 


—Bilr- 


By comparison with the exterior solution, we see that 
the interior solution has logarithmic singularities at 


[x — Bir — 1)] = 4m — 2 
and [x + Bi(r — 1)] 6B: = 4m — 2 
form = 1, 2, 3,.... These singularities are similar 


to those found by Ward! for the interior flow in a cyl- 
inder and by Pack® for the symmetric supersonic jet 
in a supersonic stream. They tend to invalidate much 
of the solution, for they appear to be a result of the 
process of linearization. 

To obtain numerical values from the above equations 
is an extremely difficult process. However, approxi- 
mations of g;’ and g»’— (a) near the jet exit and 
(b) far downstream—may be obtained from Laplace 
transform theory. 

(a) Consider the value of g.’ near x — Bo(r — 1) = 0. 
For p very large and Real p > 0, we have 


[(1 + k)/2p-V/ Bio] e? 
2TiBep Bri (ABe) 


where J’(// 8) is the function defined and tabulated by 
Ward.' Now 


= 


V 


— pBo 
we 


Ki => 
V 2 
Hence, substituting these approximations into Eqs. (15) 
and (16) yields 


B(1 + k)V/r 


[1 + O(e7 7?) ] 
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For r = | the inclination of the flow is given by 


dye? Or = [1 + (Ogo Or) + (Og2” Or)] a cos 6 = 
la cos OV (//B2)]/(1 + k) 


This can be compared with the results from two-dimen- 
sional small perburbation theory since, in a small re- 
gion near the rim, the flow may be considered two- 
dimensional. Let P’ be the pressure at the rim of the 
tube, and P? the pressure on the jet interface, then, at 
0, 


(7? Pp’) (1 2) = Bo 
(P P) (1 2)y MPP = — (2Av By) 


where Av is the inclination of the stream. Now from 


Ward's solution at the rim we have 

(P’ — Ps) (1/2)y2MeP2 = — [2a cos 0V(1/B2) /B2] 
Eliminating P’ gives 
(P Ps) (1 2) Po = — [2a cos Va B») 

Be) 
(P — Py) = —(2Av 

Eliminating P gives 

Av = cos (1 + Rk)| + [Boe (1 + R)] 
The contribution due to the angle of attack agrees 
with Oy.'*’ Or given above, and the last term agrees 
with the results of the symmetrical solution, Eq. (10), 
when v + — 1) is small andr = 1. 

(b) The asymptotic value for large x can be com- 


puted by approximating the transforms for small p. 
Hence, 


Dp) = Ti’ — (pp2) 
= (1 (1 + 61) 


We shall consider two cases, /—~ and/—~0Q. Asx 
becomes infinite, gy» becomes | 7, and we have 


dp 


r J Br: (ABe) 


j _f \ dar 
— + &) 


(1 [e"dp/ p(x — p)] = (1/d) — (e* A) 


Now 


Substituting in the above equation yields 
= (Lr) + [1 (1 + 1) — go(x + 


As/— ©, we obtain 


Ky (Ae) 1 
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go(x + i; 


Thereiore, for / infinite, the jet far downstream behaves 
like a solid cylinder at an angle of attack. If we let 


x— and/—Q(0, then we obtain 
go(l, 1) ~ 0 and = (1 + (1 r) 


Hence, the far downstream inclination of the jet whose 
diameter is large compared with the tube length is 


Av = }1 — [&/(1 + a@cosé 


Av = acos@/(1 + 


One might expect for finite length of tube that the far 
downstream jet inclination would vary between 0 and 
acos@ (1 + 


FIRST-ORDER APPROXIMATION 
It can be seen from Eq. 14 that, if we let / > 0, 


As previously, assuming small differences in pressure 
and Mach Number, we obtain, after combining g»” and 


2ri 
Substituting 


Ii(s) = nf cos 6 dé 
0 


and interchanging the order of integration yields 
¢g2’ = (1 mr) X 
x or cos '(r—x/p) 
cos V (x B+ cos 6)? — r? dé 
0 


As in the symmetrical solution, this also is the exact 
linearized solution for k = 1 and 6; = pb». 

If we consider r large (say greater than 4 or 5), we 
may write 


K,(sr) = re" V 2asr 


FIRST ORDER APPROXIMATION, LARGE 


+— -| — 
4 8 1.2 L+ te is 20 22 
| 
| 


>. 5. Slope of streamlines due to angle of attack = (1 + 
gor )a cos 6. 
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and ¢»’ becomes 
Bry 

This readily takes the form 

go! = (W2/r/x) X 
or (r—x/p) 

{ cos 6 Vx B—r-+ cos dé 

0 


or, after integrating by parts once, we obtain 
= (1/rv/2r) X 
x or cos” 
(sin? 6 Vx B — r+ cos 6) dé 
0 


= (1/rv/2r) X 
r or cos” '(r—x/B) 
Vx B — r+ cos@) — 
70 


or cos”! (r—x/B) 
(1 ror) [cos? (x B — 


r + cos x B—r+cos6d9 
Now 
cos? 6/(x8 — r + cos 6) =cos @ — (x/B — r) + 
(x — r)?/[(x/B — r) + cos 
then becomes 
(3/2)¢2’ = V2r) { [1 — (x/B — X 
or cos '(r—x/B) 
V (x B — r) + cos@+ 
0 
rorens '(r—x/B) 
(x/B — nf V (x B—r) + cos do} 
0 


For (x/8 — r) > 1, the upper limit is 7 and the in- 
tegrals are readily reduced to elliptic functions by 
substituting 


cos @ = 1 — 2sin? 6/2 = 1 — 2sin? ¢ 
and m2? = [(x/8 — r) + 1]/2. 
becomes 


= (4m (201 — + 
— 1) E(1/m)] (18) 


for m2” = [x — B(r — 1)]/28> 1. The variable upper 
limit must be treated by substituting 


The velocity potential 


= 1—2sin?¢ and sin ¢g = msint 
The velocity potential then takes the form 
(3/2)go’ = (1/wv/2r) — (x/B r)*] X 
[dt V2V1 — sin? t] + (x/8—1r) 
0 
(n2” cos? Wi - ne? sin? f)dt! 
0 


which finally reduces to 


= (4/347) [(1 — + 
(2n2* — 1)E(m)] (19) 


JOURNAL OF THE AERO/SPACE SCIENCES 


AUGUST, 1058 


where no” = [x — B(r — 1)]/28 


Eqs. (18) and (19), of course, also represent the angle 
of attack solution for /—> 0, = Mzand = 

Using the formulas on pages 73 and 76 of reference 7, 
we find that 


dK dx = |E 2x(1 — x)] — (K/2x) 
dE/dx = (E — K) 2x 


where x = for < 1. Therefore, differentiating 


with respect to 7, we obtain 
go,’ = (1 avr) [K(m) — 2E(m)| 
for 0 < m <1. Similarly, 
gor’ = (1 [(2m? — 2) 92 — 


for mn > 1. Fig. 5 shows these relations. It is noted 
that to a first approximation the deflection of the 
streamlines due to angle of attack is a function only of 
the mean Mach Number of the two flows and does not 
depend upon the pressure difference. It should also be 
pointed out that ¢g», becomes infinite along the Mach 
line x — B(r — 1) = 28, as in the symmetrical solution. 


THE SYMMETRIC Two-DIMENSIONAL FLow 


The symmetric two-dimensional flow was calculated 
by Pai,” Hawamura,' and Pack.’ For the sake of com- 
pleteness, the Laplace transforms are given here 


+ Rk) (1 — pe 


—2 — pBo(y—1 
el —e PBI) — 1) 


p(1 +k) (i — we 


where » = (1 — &) (1 + R). The potentials can be 
represented by an expansion of simple waves by first 
expanding the transforms in a series of e 7/8 and in- 
verting term by term. For any given x and y the po- 
tential will consist of only a finite number of terms 
which can be summed. The result for the interior 
solution is 


[Bre Bi(1 k)| (1 = u)] x 
(x + By) + = (x — 


where the terms multiplied by « + @,y are valid for 


x + By > (2m, + 
< (2m, + 


and the terms multiplied by (x — §,y) are valid for 


(x — By) > (2m2 + 
< (2m2 + 


where m, = 0, 1, 2, 3, etc. 
and me = 0, 1, 2, 3, etc. 


The various regions defined by the values m, and m, are 
shown in Fig. 6. Similarly, the exterior flow becomes 


Jeu” (1 + k)] (x = Boy) 
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for y > | and for 


+ 1); Bo(y — |) > 
m = GO, 1, 2, 3, ete. 


\ solution of this form was first obtained by Pai? and 
later by Pack.’ According to Pai’s theory, however, 
the values .1/; and .\/, correspond to the Mach Numbers 
where the pressures of the two streams are equal. This 
accounts for the different expression for yw. It is ap- 
parent from the form of the preceding equation that 
the slope of the boundary Og» Oy oscillates as we pro- 
gress downstream when wp < 0 and does not oscillate 


when > O. 


THE ANGLE OF ATTACK SOLUTION FOR THE Two- 
DIMENSIONAL JET 


In two-dimensional flow, the region —/ <x — 8(y — 
1) < 0 for y > 1 has a uniform velocity from the ex- 
pansion through an angle a around the leading edge of 
the upper plane. Similarly, the flow in the region 
~/ <x + @(y + 1) < 0 results from a compression at 
the leading edge of the lower plane. With this, we can 
immediately deduce the perturbation potentials de- 
finedin Eq.6. Hence, 
= — ¥ = 
l<x —B(y-—1)<0, y2 1 


Bly + 1) <0, — 


II 


If we let 
= got 


for x — B(y — 1) > O, the boundary conditions in 
Laplace transforms become, for y = +1, 


Sly P2y 
For 2; and @.” we assume 


= fi(p) sinh 


"=g(ple for y> 


Gi 
te 
II 


= go(pre' for ys 


Substituting into the boundary conditions, we find 
first that go = —g,; and after some calculation, we 
finally obtain for the transforms 


sinh pp.y 
p*Bi(cosh + k sinh 


¥1 


cosh pp; 


(cosh + k sinh pps)’ 


Replacing the hyperbolic functions by exponentials 
and expanding in simple waves yields the following ex- 
pression for gy’: 


Fic. 6. Two-dimensional flow. Model and coordinate system 


my 


p + — 
+ tu 0 


x + Ay) — 


which is valid for x and y satisfying 


x + By > (2m + 
< (2m + 3)61 


x — By > + 
< (2m. + 


Similarly for ¢g.”, we have 


BL R) l+u 
1—(- ym 


which is valid for y > 1 and 
x — Bly — 1) > 2mp, 
< 2(m + 1)6; 
Combining with ¢g) and simplifying by noting that 
1/(1 + k) = (wu + 1) 2, we have 
= (1 — pw) [x — Bly — 1)] 
Similarly, for y < — 1, 
= (1 — pw) + + 1)) 22 
which is valid for 


x + + 1) > 


< + 1)8,, = 0, 1,2,3,... 


The inclination of the jet boundary is found to be 


all + (Oge’, Ov)] = all + (u — 1) 


which is seen to depend only on wu and a. If uw = 0, 
then the flow is deflected only in the interval 0 < x — 
Bo(y — 1) < 28, and becomes parallel to the stream for 
x > 26;. The second term always has the same sign 
if uw is negative —that is, ifk > 1. Fork < 1, the sign 
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of the second term alternates with consecutive values of 
m, and the boundary has a waviness whose amplitude 
approaches zero as we move downstream. We also note 
that when the symmetrical solution shows waviness, 
the angle of attack shows none; while if the slope of the 
symmetrical solution is monotonic, the angle of attack 
term shows waviness. 


CONCLUSION 


Solutions of the linearized equations for the flow of 
a supersonic jet at an angle of attack in a supersonic 
stream have been found for both three-dimensional and 
two-dimensional flows. The symmetrical three-dimen- 
sional jet has a boundary which has an oscillating slope 
as we move downstream from the exhaust when k = | 
where 


V Me? — 1 


k= 
¥2P VME 1 


The frequencies are given by the imaginary part of the 
complex roots of the equation 


D,(p) Ko( pee) Ki = 0 


The perturbation due to the angle of attack also shows 
a waviness. For this part of the solution the fre- 
quencies are given by the imaginary parts of the com- 
plex roots of 


= — pa.) = 0 


For k = 1, yifi = y2p2, and MJ, = Mz, the jet interface 
has no periodicity, but the downstream Mach cone 
from the axis at a point a distance 8 = VM: — 1 
from the jet exhaust becomes a surface of discontinui- 
ties where the velocity components become infinite. 
In general, the potentials for both the symmetric and 
angle of attack flows exhibit unrealistic singularities 
along certain Mach lines in the flow when k + 1. 

The length of tube before the cylindrical jet exhaust 
has some influence on the deflection of the jet due to the 
angle of attack. When the length of the tube is in- 
finite, then far downstream the jet behaves like a solid 
cylinder. If the length of the tube is zero, then the 
inclination angle of the jet far downstream is 


a [1 + y2P2M,2")] 


where a is the angle of attack of the outside stream. 
The inclination of the two-dimensional jet far down- 
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stream is equal to a and does not depend on the length 
of the tube. 

For the two-dimensional jet the parameter 2 is a dj. 
rect measure of the reflection of waves from the jet 
boundary. For k = 1, no reflections occur and the 
symmetrical jet expands to a uniform flow in a short 
distance of 28,. Similarly, for the perturbation due 
to the angle of attack, the flow aligns itself with the 
stream in the distance 23;. The role of the parameter 
k in the reflections irom the boundary was noted by 
Wilder and Hindersinn,‘ who stated the parameter in- 
correctly in the form 


pil p2l 


where p is the density of the stream. It should have 
been written 


pil pol 


which is the exact equivalent of our definition for k. 

For k > | and zero angle of attack, the boundary of 
the two-dimensional jet oscillates between positive and 
negative slope as we move downstream, while, for k < 1, 
the boundary slope decreases monotonically. For both 
cases, the magnitude of the slope decays exponentially, 
The angle of attack term, however, indicates a waviness 
for k < 1 and a monotinically decreasing slope for k > 1. 
The wavelength of the oscillations for the two-dimen- 
sional jet does not depend upon the parameter & but 
is equal to 49). 
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Buckling Phenomena of Stiffened Panels’ 


SYED YUSUFF* 
Bristol Ldircraft Limited 


SUMMARY 


An investigation of the instabilities involved in the design of 
It is 
shown that three modes which constitute predominantly a skin 


stiffened panels forming compression surfaces is presented 
buckling phenomena are largely encountered in compression 
panels—the first mode, which is associated with a long wave- 
length (same wavelength as Euler instability), is due chiefly to 
flexural deformation of the stiffeners; the second mode, which 
occurs in short wavelength (length of buckle is of the order ot 
stiffener spacing), results mainly from the torsional restraining 
effects of the stiffeners; and the third mode is a combination of the 
first two modes with the latter of these occurring as an initial 
instability and the former at failure. 

The theory presented, although general, is well illustrated 
by the two most popular and commonly used methods of con- 
struction—-namely, Z-stiffened panels of riveted construction and 
integral construction with unflanged stiffeners. 

An extensive analysis of test panels is carried out; these panels 
cover a wide range of configurations required to satisfy wing 
loadings typical of modern aircraft. The results of this analysis 
are found to be in good agreement with the test data. 

Tabulated results of the analysis and charts to determine criti- 
cal stresses are offered as a guide to design. 


SYMBOLS 
A, = stiffener area 
b = stiffener spacing 
Cc = torsional rigidity of stiffener 
D = flexural rigidity of plate of unit width = 
| Et8/12 (1 — v?){ 
ds = attached flange width 
dr = outstanding flange width 
E = Young’ modulus 
G = shear modulus 
h = height of stiffener 
f, = effective moment of inertia of stiffener 
L = pin-ended length of panel 
m = number of half-waves 
t = thickness of skin 
: = thickness of stiffener 
w = deflection of skin or stiffener 
9 = coordinate axes 
= compression stress 
Oe; = critical compression stress 


p = radius of gyration of the section consisting of 
a stiffener attached to the center of skin of 
width } 

\ = length of buckle (= L/m) 

Poisson's ratio 

nondimensional parameters defined in the text 


plasticity correction factor 


INTRODUCTION 


I ORDER TO WITHSTAND high loadings and at the 
same time to maintain good aerodynamic smooth- 
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ness, stiffened panels forming the compression surface 
of a wing are designed to high initial instability stress. 
Furthermore, in order that the surface smoothness is 
undisturbed, it is desirable to proportion elements of the 
panel in such a way that initial instability and failure are 
coincident. Such panels should also be designed for 
minimum weight. 

These requirements can be satisfied by panels in- 
volving either thick-skin construction or closely spaced 
stiffeners (15 < 6 ¢ < 35 for high wing loading), which 
may be riveted or integrally machined with the skin. 

It is well known that the solution of stiffened panels 
meeting all these requirements is rather complex. 

Various authors'~* have studied the instability of 
open section stiffeners of stiffened panels under uniform 
compressive stress. Their investigations have been 
essentially the extension or adaptation of torsion-bend- 
ing theory, ® the effects of the sheet in restraining 
twisting and warping of the stiffener as well as in fore- 
ing the center of the twist to lie in the plane of the sheet 
having been considered. This instability was con- 
sidered to consist of twisting of the open-section stif- 
feners about an axis parallel to their centroidal axis 
without cross-sectional distortion. (The type of fail- 
ure goes from combined twisting and bending to pure 
bending as the axis of twist is moved from the neigh- 
borhood of the stiffener to infinity.) 

Levy and Kroll* have carried out a somewhat ex- 
tensive comparison of their method with experiments. 
For panels reinforced by Z-stiffeners and having a 
stiffener thickness of two or more times the sheet thick- 
ness, they observed that the agreement was good; 
however, in panels having stiffener thickness equal to, 
or less than, the sheet thickness, their proposed method 
of computing stiffener loads yielded total loads up to 
35 per cent too large. 

Because of such large discrepancies, it is likely that 
torsion-bending mode may not be present in panel hav- 
ing the sheet thickness equal to, or greater than, the 
stiffener thickness, and that failure may be due to 
other buckling phenomenon. 

Hence, in the present paper, a different approach, 
which was first considered by the writer’ ~* for the study 
of integral construction with unflanged stiffeners, is 
used to investigate the buckling phenomena of stiffened 
panels. It is found that possible buckling modes, 
apart from that of Euler, which may be responsible for 
failure of panels, are chiefly due to (a) cross-sectional 
distortion of the stiffener in its own plane and asso- 
ciated interaction with the sheet while the sheet- 
stiffener junction remains straight, and (b) flexural de- 
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Fic. 1. Z-stiffened panels; symbols for dimensions. 


Fic. 2. Integral construction; symbols for dimensions. 


formation and rotation of the stiffener without distor- 
tion in its own plane of the cross section. 


The instability associated with the deformation 0° 


type (a) is known as local (or crippling) instability. 
For integral construction with unilanged stiffeners, the 
local instability stress can be determined using the 
chart given in reference 10; but in cases of panels rein- 
forced by Z-section stiffeners, local instability of 
stiffeners is not of much significance in almost the en- 
tire range of practical stiffened panel design (see Tables 
1 to 3). Therefore, the references to various studies 
dealing with this instability in Z-stiffeners are omitted. 

The analysis presented here deals mainly with the 
instabilities associated with deformations of type (b). 
Two primary modes of instabilities arise, depending 
upon whether the restraints offered by the stiffener on 
the sheet are predominantly flexural or torsional. Be- 
cause of the wavelengths associated with these types, 
they are referred to herein as long-wave buckling and 
short-wave buckling modes. In both these cases, a 
plate buckling phenomenon is involved. It is found 
that the long-wave mode is of great interest and, in 
fact, constitutes the basis for most designs. 

Although the analysis presented here is quite gen- 
eral, particular reference is made to Z-section stitened 
panels (Fig. 1), and to integral construction with un- 
flanged stiffeners (Fig. 2). Because of their high effi- 
ciencies and, also, in the latter case, of fewer machining 
difficulties, these two types are commonly used. There- 
fore, an extensive comparison of the theory with tests 
covering these two types is made. 


ANALYSIS 


In the development of the analysis the following as- 
sumptions are made. 

(1) The stiffeners are regarded to have sturdy cross 
sections; that is, stiffeners which do not suffer cross- 
sectional distortions in their own planes. This assump- 
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tion is reasonable for the type of Z-stifener panels 
used in aircraft construction. However, in the case 
of unflanged stiffeners of integral construction, the 
height of the stiffener is fixed so as to avoid local in- 
stability. 

(2) Each stiffener, having flexural rigidity //, and 
torsional rigidity C, is assumed to be concentrated at 
the sheet-stiifener junction in the plane of the sheet. 
Its effective moment of inertia /,, and torsional rigidity 
C, will be defined later. 

(3) The stiffeners are simply supported at the ends 
and they are compressed together with the plate and 
-arry the same stress as the plate at buckling. 

(4) During buckling the stiffeners deflect laterally in 
the same direction and rotate alternately as shown in 
Fig. 3. The deflected surface of the plate is symmet- 
rical with respect to the x axis. If we assume that the 
stiffeners are deflected during buckling in opposite 
directions, we obtain a transcendental equation which 
always gives larger values for critical stresses than those 
obtained from the assumption of symmetry. 

(5) The material is homogeneous and isotropic and 
obeys Hooke’s law of deformation. 

Consider the sheet between two consecutive stiffen- 
ers (Fig. 3). Because of rigid connection with the 
stiffeners, it is assumed to be elastically supported with 
elastically built-in edges at the junctions. The lateral 
deflection of the sheet at the junction is not zero but 
equal to the deflection of the stiffeners. 

The differential equation of the deflected sur‘ace 
for uniform compression along the x axis is given by'! 
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Fic. 3. The panel with buckled sheet. 
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(Q'w + + (O'w Oy*) = Substituting Eq. (2) in Eq. (1), we have the diifer- 
ase —(¢,t/D) (O°w/Ox*) (1) ential equation 
the Assuming that under the action of compressive forces — + 
the plate buckles in m sinusoidal half-waves, the solu- (o,t D) (m2x? L2)]f = 0 (3) 
tion of Eq. (1) may be written Fe 
nd where 
at w = f(y) sin (mmx L) (2) as 
et f’ = Of(y)/Ox, f” = O*f(y) Ox’, ete. . 
ies in which /(y’) is a function of y alone and is to be deter- 
7 mined later. Since, owing to constraints along the sides y = +0 2, 
Eq. (2) satisfies the boundary conditions along the o,t D > m*x* L? is always true, the general solution of 
ds ends of the panel (x = 0,x = L). Eq. (3) can be represented in the following form: - 
nd 
in TABLE 1 
in Long Wave Buckling 
Ay 
Fail- 
he Crit. ing 
te Serial Length* Stress Stress 
oh Num- (2L) ts (Theory) (Test) 
ber (in. ) (in. ) t/t b/t h/t dr/t, da/ts L/p (in.*) EI. /bD (ksi) (ksi, 
” Material 75 S-T6 (Test Ref.: panels 1-6, ref. 16; panels 8-31, ref. 15) 
] 41.2 0.0637 0.252 14.9 20.1 7.89 20.37 62.5 0.1072 19.2 2.40 23.7 23.1 
id 2 28.2 0.0637 0.251 14.8 19.9 8.42 20.19 42.5 0. 1067 18.9 3.42 49.0 43.7 
3 53.0 0.0655 0.261 29.9 29.5 11.65 19.97 62.5 0.358 33.0 4.70 22.6 22.8 
4 36.2 0. 0665 0.273 30.9 29.1 11.55 19.52 42.5 0.364 36.4 6.10 35.7 32.9 
- 5 23.4 0.0663 0.271 30.7 29.0 11.67 19.74 27.5 0.358 35.5 6.30 39.0 39.6 
6 68.8 0.0637 0.261 41.0 40.2 15.91 20.38 62.5 0.767 51.0 5.73 17.9 17:8 
7 39.5 0.1018 0.410 20.2 7.95 12.7 62.5 0.564 81.0 3.1 21.6 24.8 
h | & 26.8 0.1015 0.407 20.1 20.1 8.06 12.91 42.5 0.559 80.0 4.5 165.0 43.5 
| 9 73.4 0.100 0.40 25.0 20.4 8.16 12.86 62.5 0.583 65.1 3.83 21,5 21.8 
al | 10 19.8 0.0987 0.396 25.0 20.7 8.17 13.03 42.5 0.563 64.5 5.40 12.8 46.0 
it 11 70.0 0. 0988 0.401 30.1 21.0 8.20 13.11 62.5 0.618 61.0 4.6 21.2 20.8 
12 17.0 0.0985 0.392 30.3 21.0 8.10 18.10 42.5 0.611 55.6 6.16 38.0 44.8 
13 78.0 0.101 0.409 30.2 30.1 11.92 12.74 42.5 1.718 166.5 6.48 42.0 39.6 
‘e 14 66.0 0.0997 0.654 20.7 20.6 8.11 10.01 42.5 0.433 422.0 4.58 44.8 46.7 
15 93.5 0.0998 0.65 25.4 20.4 7.96 9.78 62.5 0.481 372.0 3.95 22.2 23.7 
16 63.0 0.0993 0.642 25.3 20.6 8.07 9.98 42.5 0.426 324.0 5.4 41.8 44.7 
| 17 148.5 0.1005 0.646 25.0 30.7 12.03 9.71 62.5 1.414 1006.0 3.9 22.3 21.6 
18 101.0 0.1004 0.646 25.0 30.7 12.04 9.87 42.5 1.215 906.0 5.28 410.7 38.8 
19 47.9 0.0997 0.642 30.2 12.3 4.77 10.48 62.5 0.13 81.2 4.67 22.0 24.8 
~ 20 32.4 0.10388 0.673 30.2 1.7 4.59 9.45 42.5 0.12 77.2 6.24 39.2 43.1 
21 89.4 0.0991 0.641 30.0 20.5 8.06 9.95 62.5 0.494 317.0 4.66 22. ¢ 25.9 
22 60.6 0.1004 0.646 30.1 20.2 8.05 9.77 42.5 0.45 284.0 6.45 42.0) 47.4 
23 145.0 0.1014 0.649 30.2 30.4 11.84 9.72 62.5 1.494 898.0 4.65 21.8 23.2 
24 98.2 0.1004 0.654 31.4 30.6 11.99 9.87 42.5 1.29 806.0 6.6 40.5 39.9 
25 83.0 0. 0997 0.651 40.9 20.6 8.32 9.91 62.5 0.5415 267.5 6.4 22.6 24.3 
26 137.0 0.1012 0.664 40.8 30.4 11.91 9.71 62.5 1.615 800.0 6.48 22.8 23.9 
27 ef 0.102 0.979 20.1 20.1 7.9 6.63 42.5 0.301 1390.0 4,21 40.2 37.7 
28 106.0 0.1018 0.988 20.0 19.9 7.92 6.64 62.5 0.3445 1670.0 3.05 21.3 20.6 : 
29 71.4 0.0997 0.981 25.0 20.8 7.94 6.63 42.5 0.32 1315.0 5.27 40.7 41.5 : 
30) 104.5 0. 1008 0.982 25.0 20.4 8.0 6.51 62.5 0.372 1475.0 3.85 21.7 22.8 
31 69.8 0.0993 0.96 29.6 20.5 7.9 6.71 42.5 0.335 1090.0 6.25 41.0 42.0 
Material 24 S-T3 (Test Ref. 17) { 
32 72.6 0.064 0.51 35 24 6.2 re 76.6 0.1405 177 4.15 13.5 13. 
33 51.2 0. 064 0.51 35 24 6.2 7.4 54.0 0.1405 177 5.78 26.2 26.9 
34 60.3 0.064 0.68 35 19 7.5 10.4 76.0 0.0739 216.5 4.5 15.9 15.4 
a) 11.0 0. 064 0.63 35 19 7.5 10.4 51.75 0.0673 197 6.4 29.0) 28.3 a 
36 94.6 0.064 0.68 35 29 11.5 10.4 71.5 0.225 661.5 4.8 18.1 18.5 ae 
37 67.5 0.064 0.63 35 29 11.5 10.4 51.0 0.20 587 6.37 29.0 28.1 
38 64.0 0.064 0.79 35 19 7.5 9.3 73.0 0.06525 4172.5 4.72 17.5 16.8 : 
39 16.4 0.064 0.79 35 19 7.5 9.3 52.75 0.0579 420.0 6.08 29.1 29.7 ‘ 
40) 103.5 0. 064 0.79 35 29 6.2 5.0 72.0 0.1968 1420 4.7 17.4 17.5 : 
41 75.0 0.054 0.79 35 29 6.2 5.0 52.0 0.1688 1220 6.03 28.4 28.2 
42 85.0 0. O64 0.79 35 24 6.2 6.03 73.5 0.1203 871 4.65 17.0 7.3 
48 63.0 0. 064 0.79 35 24 6.2 6.03 54.5 0.105 760 5.85 27.0 29.3 
44 70.64 0.054 1.0 35 19 9.5 8.1 72.0 0.0624 1160 4.75 17.7 17.5 
45 50.6 0.064 1.0 35 19 9.5 8.1 51.5 0.0536 992 6.08 28.0 27.3 
Integral Construction with unflanged stiffener, Material DTD. 683 (Test Ref. 8, 18, 19) ad ei 
35.5 0.141 0.141 25.0 10.3 38.8 0. 0S09 352 6.10 58.0 58.2 
47 51.6 0.183 1.785 29.4 9.95 43.8 0.1905 641 6.38 46.0 18.6 
48 51.6 0.2017 1.97 29.4 9.3 : ; 42.0 0.2163 740 6.59 50.0 52.0 
0. 1432 1.40 33.7 14.0 = ‘ 47.0 0.218 643 6.83 40.0 33.1} 
0.1451 1.465 35.3 12.1 ; : 48.2 0.1458 490.5 6.80 36.8 37.4 


* Ends assumed to be clamped; the fixity coefficient for tests at NACA = 3.75. 
: +t E = 10.91 X 10° psi, vield stress 73.7 ksi; limit of proportionality 50.0 ksi (average values determined on specimens extracted from 
six slabs ) 

t Panel failed by local instability (see reference 7). 
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Serial 
Number 
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
9 
26 
27 
Serial Length* 
Num- (2L) 
ber (in.) 
1 30 
‘ 2 71.5 
3 56.2 
4 67.6 
5 46.9 
6 48.0 
i 7 80.0 
; 8 109 
9 60.2 
10 95.5 
11 130 
12 40.7 
13 128 
14 57.8 
15 121.5 
16 90.0 
17 109 


* Ends assumed to be clamped; fixity coefficient 


ty) 


where 


The constants of integration in this solution can be 
determined from the conditions of the constraints along 
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Short Wave Buckling, Material 75 S-T6 


SPACE SCIENCES—AUGUST, 1958 


the sides y = +) 2. For the case under consider. 
ation they are obtained by the requirements that the 
rotation of the edges of the sheet is equal to the twist- 
ing and their deflection is equal to the deflection of the 
stiffeners. 

To obtain the first boundary condition the twisting of 
the stiffener is considered. The angle of rotation of 
its cross section is Ow Oy, and the rate of change of 
this angle along the edge is O°w/Oxdy. Hence, 


TABLE 2 


Critical Stress 


Length* —— ~ 
(2L) ts Theory Test 
(in. ) (in.) b/t h/ts dr/ts d4/ts (ksi) (ksi) 
(Test Ref. 16) 
21.1 0.0665 0.276 41.6 29.1 11.55 19.67 20.9 20.8 
32.2 0.0673 0.273 40.6 28.5 11.38 19.29 21.9 21.5 
46.9 0.0656 0.270 41.2 39.0 15.45 19.79 21.3 22.5 
63.0 0.0644 0.264 §1.1 39.8 15.81 20.23 13.9 14.2 
39.2 0.0658 0.267 61.0 29.6 11.68 19.81 9.7 9.9 
26.5 0.0656 1.028 50.6 29.2 ‘1 74 6.6 20.3 21.8 
42.0 0.0658 1.045 51.0 29.1 11.67 6.5 20.2 21.4 
36.4 0.0644 1.014 50.6 39.6 15.80 6.49 20.8 19.3 
57.0 0.06438 1.007 50.0 39.9 16.17 ef 21.3 20.3 
(Test Ref. 15) 
32.2 0.0982 0.394 25.0 20.9 8.21 13.20 57.0 56.7 
20.0 0.1018 0.406 25.0 20.0 7.96 12.99 57.0 8 4 
30.2 0.1018 0.41 30.0 20.0 7.89 12.72 40.2 43.0 
18.8 0.102 0.412 30.1 20.0 7.9 12.92 40.0 41.4 
16.25 0.0986 0.406 30.1 20.1 8.11 13.18 40.0 40.6 
14.8 0.1016 0.653 25.0 20.2 7.93 9.7 63.0 61.2 
26.15 0.0980 0.634 25.2 20.8 8.12 10.01 63.0 61.2 
14.2 0.1018 0.656 30.3 20.0 7.94 9.58 45.0 48.4 
24.7 0.0991 0.635 30.2 20.5 8.08 9.7 45.0 47.4 
23.7 0.1009 0.656 30.0 30.6 11.95 9.77 48.0 43.6 
40.4 0.1016 0.654 30.1 30.2 11.86 9.65 47.8 44.9 
21.7 0.10038 0.657 41.0 30.5 12.01 9.72 24.8 24.9 
38.1 0.1021 0.666 40.9 30.1 11.82 9.62 25.0 2.1 
9.8 0.0977 0.957 29.5 12.7 4.9 6.91 55.7 57 3 
16.4 0.0982 0.946 29.1 12.6 4.9 6.80 56.0 54.0 
16.8 0.1024 0.998 29.6 20.0 7.87 6.78 58.0 54.6 
28.8 0.10038 0.973 29.4 20.4 7.96 6.64 58.0 56.4 
16.0 0.1017 0.981 39.3 20.2 7.91 6.65 32.4 33.8 
27.6 0.0995 0.963 39.3 20.6 8.1 6.90 32.0 32.5 
TABLE 3 
Double or Mixed Instability, Material 75 S-T6 
Short- 
Crit. ening 
Crit. Stress stress Av. per unit 
(shortwave) (long- failing length 
v - wave) stress at fail- 
ts Ie (long Theory Test (theory) (test) ing load 
(in.) b/t h/ts d4/ts L/p El,/bD wave) (ksi) (ksi) (ksi) (esi) 
(Test Ref. 15) 
0.1018 0.41 30.0 20.0 7.89 12.72 27.5 0.605 58.3 6.82 40.0 413.0 47.2 49.3 510 
0.101 0.409 40.5 30.4 11.91 12.82 42.5 1.893 125.5 7.62 22.0 24.4 32.4 34.2 330 
0.0991 0.652 41.1 20.6 8.11 9.91 42.5 0.488 244 8.33 24.2 27.6 37.7 39.4 $30 
0.0999 0.966 39.2 20.4 7.94 6.87 42.5 0.361 892 8.12 32.4 33.0 39.3 40.3 106 
(Test Ref. 16) 
0.0656 0.27 41.2 39.0 15.45 19.79 42.5 0.798 60.6 6.95 21.3 22.5 26.0 27 266 
0.0641 0.622 19.6 19.8 7.93 9.57 62.5 0.0901 177 7.3 15.9 17.4 19.8 21 223 
0.0647 0.63 49.6 29.7 11.84 9.55 62.5 0.2769 542 7.7 16.5 17.4 22.0 22.9 230 
0.0639 0.609 58.3 40.0 15.9 9.95 62.5 0.6065 935 8.8 11.9 13.0 20.8 18.2 198 
0.0635 1.007 51.0 20.3 8.13 6.90 62.5 0.0718 975 8.0 18.5 21.4 22.6 23.2 226 
0.0661 1.037 50.3 29.0 11.62 6.32 62.5 0.211 2340 7.65 20.6 71.2 21.2 22.1 202 
0.0658 1.035 50.3 38.9 15.61 6.50 62.5 0.4315 5180 7.3 21.4 19.2 19.3 19.6 186 
0.0637 1.018 51.5 19.9 8.16 6.9 42.5 0.06 835 10.27 18.5 21.0 36.3 35.8 120 
0.0648 0.939 56.6 39.6 15.7 6.6 62.5 0.4665 3970 8.5 15.8 15.9 20.6 19.4 174 
0.064 0.967 57.9 19.9 7.86 6.45 62.5 0.0743 740 8.8 13.8 16.3 71.1 20.1 197 
0.0602 0.868 69.1 42.6 17.02 6.95 62.5 0.41 2790 9.85 9.75 10.3 18.5 17.8 183 
0.066 1.055 76.8 28.9 11.65 6.82 62.5 0.2295 2135 11.5 8.65 9.4 20.6 19.6 228 
0.0628 1.022 76.7 40.4 16.05 6.56 54.5 0.404 3820 12.14 8.9 8.8 23.0 23.3 239 
= 3.75 
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BUCKLING 


the twisting moment at any cross section of the stiffener 
in the direction of the x axis is C(O*w/Oxdy). Since the 
plate connected rigidly to the stiffener transmits con- 
tinuously distributed twisting moments, this moment 
varies along the edges of the stiffener. The magnitude 
of these applied moments per unit length is numerically 
equal to the bending moments transmitted by the plate 
on either side of the stiffener. It can be shown that 
for symmetrical buckling the moment transmitted by 
the plate on each side of the stiffener has the same 
sense. Therefore, equating the bending moment of 
the plates with applied moment per unit length of the 
stiffener, we have 


2D[(d°« Oy?) + v(0?w Ox") | = 
Ox*0y) for y= (Sa) 


—2D|(0°w + Ox*)] = 
C(0*w Ox*d0y) for y = —b/2 (5b) 


To obtain the second set of boundary conditions, 
the bending of the stiffener is considered. Since the 
pressure that is transmitted from the sheet on the stiff- 
ener on either side has the same sense, the intensity 
of load due to this is 
q = 2D|(0*w Oy’) + (2 — v) /Ox*dy) ] 

for y=b2 (6a) 


g = —2D[(0*w/dy*) + (2 — v) Ox*dy) | 
for y = —b/2 (6b) 
The differential equation of the deflection curve of 
stiffeners is given by 
EI,(0'w Ox*) = + (2 — v) X 
(O*w Ox*dy)] — A,o,(O°w Ox?) for y= 2 


EI. (O'w Ox!) = —2D[(O*w Ov*) + (2 — v) X 
(O*w Ox*dyv)] — A,o,0°w Ox? for y = —b 2 (6d) 


These two equations form another pair of boundary 
conditions. 

Using the boundary conditions (5a), (5b), (6c), and 
(6d), the four constants of integration in Eq. (4) can be 
determined. However, for the symmetric case, as is 
envisaged in the present investigation, Eq. (4) can be 
written: 


f(y) = A cosh ay + B cos By (7) 
and with this, Eq. (2) becomes 
IV = (A cosh ay + B cos By) sin (rx X) (8) 


Substituting the boundary conditions (5a) and (6a) 
in Eq. (8), we have, respectively, the following: 


A}2De? cosh a(b/2) — 2vD(m )? cosh a(b 2) + 


C(r \)? a sinh a(b 2)} —B}2D~B? cos B(b/2) + 
2vD(r cos B(b/2) + Clam d)28 sin B(b/2)} = 0 (9) 


X)4 cosh a(b 2) — sinh a(b/2) + 


2(2 — sinh a(b 2) — Ayo,(m/d)? X 
cosh a(b 2)} + cos B(b/2) — 
sin B(b/2) — 2(2 — X 


sin B(b/2) — A,o,(m dX)? cos B(b 2)! (10) 
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To determine the critical stress, the determinant of 
the coefficients A and B in Eqs. (9) and (10) is equated 


to zero. The resulting equation is given by 


[(a? + 8?) 2] [(EI, D) (wd)? — (Ayo, + (1 4) X 
(C D) (wd)? D) \)? —(A,o, D)! X 
}a@ tanh a(b 2) + B tan B(b 2)} — 
(1/2) (C/D) (a? + B?) } ag tanh a(b 2) X 
tan 6(b/2)} — a tanh a(b 2) | — 
(2 — v)B? — o(2 — v) (w dA)? va} 
tan B(b 2) [a2B? + (2 — + 

v(2 — v) — =O (11) 


To simplify this equation further, we introduce the 
following notations: 


mrb L » =bDC, D=y 
(EI, /bD) — (A, bt) o)? = 6 
Substituting these in Eq. (11), we get* 


+ + ¢) X 
tanh (1 2) V oly +o) + V oly — o) X 
tan (1 2) Vo(y — —(¥o + 6) X 
tanh (1 2) + oly — X 
tan (1 2) — — — 9) + 
(1 — vp}? tan (12) Voly — — X 
ly — (1 — vo}? tanh (1 2) +o) = 0 (12) 


This equation describes the behavior of stiffened 
panels when the interaction between the sheet and 
stiffener is such that both suffer the same rotation and 
lateral deflection at their junction. For the purpose of 
design this equation appears to be lengthy and rather 
cumbersome to solve. However, further simplification 
gives rise to three particular cases, which, without much 
loss of accuracy, mainly represent the behavior of stiff- 
ened panels. The three cases are discussed below. 


LONG WAVE BUCKLING 


When stiffened panels buckle under uniform com- 
pressive stress, the number of half-waves formed in the 
pin-ended length LZ of panels may vary from one to m. 
We are now interested in those panels that buckle so 
that there is only one-half wavelength in the effective 
pin-ended length of the panel (one wavelength in the 
case of panels loaded between flat ends or with ends 
clamped). In this case, the major deformation of the 
stiffeners is flexural. Whenever buckling of this 
type occurs, the failure of the panels is expected to be 
coincident with the onset of buckling. Most of the 
panels which fall within the range of practical design 
(Table 1) buckle in this manner, and failure takes place 
suddenly and without any warning. 

In panels which are designed for high wing loading, 
and to achieve high stress over as large a rib-pitch (un- 
supported length) as possible, it is assumed that the 


* A similar equation for the instability of a plate supported on 
two beams was first obtained by Chwalla.'* 
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Fic. 4+. Relation between buckling and aspect ratio parameters. 
flexure of stiffeners is more significant and the wave- 
length is the same as that of the Euler instability. In 
such conditions the sum of the terms in Eq. (12) relat- 
ing to the torsional rigidity is small when compared to 
the sum of other terms. Hence, by omitting these 
terms, Eq. (12) is simplified to 


— Vw — + (1 — tan (1/2) X 
Voy — — X 


tanh (1/2) + = 0 


(13) 

The solution of this equation* is represented graphi- 
cally in Fig. 4. By extending this Figure to the right, 
it can be shown that for @ = ©, the minimum value of 
the buckling parameter y is 6.28. The stress corre- 
sponding to y = 6.28, for vy = 0.3, is oo = 3.62 E(t db), 
which is the same as the critical stress of an infinitely 
long plate having its unloaded edges simply supported. 

When the torsional restraining effect of the stiffener 
on the sheet is taken into account, the minimum value 
of the buckling parameter, y, can vary from 6.28 to 
8.3 (corresponding to buckling of an infinitely large 
plate with its unloaded edges simply supported or 
fully fixed). Hence, in analysis of long-wave buckling 
which causes failure as an initial instability, the mini- 
mum value of y may be in this range (see values of y 
in Table 1). However, experiments’ on Z-stiffener 
panels having b/¢ = 30, h/t, = 20, and t,/t ~ 1 have 
shown that the maximum value of y achieved is 6.28 
over a large range of slenderness ratio (34 < Lp < 42). 

For long-wave buckling, the solution of Eq. (13) con- 
sists of finding y for any given value of ¢ (= 7d/L), 


* Compare with Eq. (t), on page 347, reference 11. 
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m = 1, by trial and error. The chart given in Fig, 4 
can be used or numerical methods employed to deter. 
mine correctly the value of ¥. The critical stress js 
then obtained from the following equation: 


= [E1201 — (Ida 
in the elastic range, and 


gen = (E11 = 


above the limit of proportionality. In this equation » 
is a plasticity correction factor. In plastic buckling of 
panels of integral construction, critical stresses have 
been computed from Eq. (l4b) by using = 
(Ey = tangent modulus). The results are consistently 
found to be in good agreement with the tests. 


EFFECTIVE MOMENT OF INERTIA 


In Eq. (13), the parameter 6 contains /,, the effective 
moment of inertia of the stiffener. This must be known 
in order to solve this equation. Because of the large 
number of sizes and shapes that are commonly used in 
both riveted and integral construction, the interaction 
between the skin and stiffener is naturally complex and 
determination of the effective moment of inertia is not 
readily obvious. It was suggested by Timoshenko,'! 
for problems in stability of plates reinforced by ribs in- 
volving large width of plate and stiffeners such as Z 
and channel sections, that the moment of inertia of! the 
stiiTener should be calculated about the axis at the inter- 
face, parallel to the plane of the sheet, and, due to the 
rib being rigidly connected with the plate, a portion of 
the plate should be included in obtaining the flexural 
rigidity of the rib. 

On the basis of similar considerations, test results, 
and the nature of interaction between the skin and 
stiffener, the effective moment of inertiaT of the stiffener 
is defined as follows. 

For computing the effective moment of inertia, an 
element of the panel consisting of a stiffener at the 
center of the skin corresponding to the stiffener spacing 
is considered. The effective moment of inertia, /,, for 
the two types of construction considered here is given 
by (a) integral construction with unflanged stiffener 


3[1 + (1/2) 
bt \ [1 + (A, 


obtained by adding the moment of inertia of one-half 
of the skin area |(1 2)d¢] about the centroid of the ele- 
ment to the moment of inertia of the stiffener about the 
same point, and (b) Z-stiffened panels [(b ft) < (h ,)], 
for panels having ¢, ¢ > 0.6 and 40 < L p < 62.5, and 
t,/t>0.6 and 54 < L/p< 40, 


I, = 1, + {A, ¥,2/[1 + (A,/bt)]}?} (16) 


in which /, is the moment of inertia of the stiffener about 
its centroid and Y, is the distance between the centroid 
of the stiffener and the plate middle surface. It may be 


+ For further information about the effective moment of inertia 
of a stiffener, see reference 13. 
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BUCKLING 


noted that the effective moment of inertia given in Eq. 
(16) is the moment of inertia of the stiffener about the 
centroid of the plate-stiffener element. 

For panels having f, f > 0.6 and L p= 60, 


1A, {1 + 
(1/2) 


[1 + (A, bt)]*} (17) 


This equation is obtained by adding the moment o: 
inertia 02 one-half of the skin area [(1 2)/t| about the 
centroid of the plate-stiffener element to Eq. (16). 

It may be noted here that the dimension of the out- 
standing flange, d,, is fixed from efficiency consideration 
and is 0.44, whereas the width of the other flange, d4, 
is determined by a desirability of suitable attachment 
such as rivet diameter, rivet pitch, etc. 


SHORT WAVE BUCKLING 


This refers to buckling of the plate owing to the 
twisting of stiffeners. 

In panels having small slenderness ratio, or large 
values of 6, the sheet-stiffener junction is assumed to be 
straight. Thus, there are no flexural deformations of 
the stiffeners, and only their torsional restraining effects 
need be considered. Buckling in the short wave mode 
can therefore occur earlier than that in the long wave 
fashion. Since in the short wave mode, in which half 
a wavelength is of the order of stiffener spacing, ¢ ~ 7, 
it can be shown that the first two terms in Eq. (12) are 
very large when compared to the remaining terms. Re- 
taining only the first two terms in the equation, we have 


+ (1 4)640/n) + X 
tanh (1/2) Vo(y + 6) + — X 
tan (1/2) Vd(y — ¢)} =O (18) 


This equation* is simplified further to 


/ 
+ ¢) tanh (1 2) Vo(y + 6) + 
V o(y — ¢) tan (12) — ¢) + 
+ (uy = 0 
The graphical solution of this equation is given in Fig. 5, 
which is drawn from charts given in reference 14. The 
nondimensional parameter yu, for various values of which 


the chart in Fig. 5 is drawn, is defined as 


p = OD/C = bD/GIJ 


(19) 


(20) 
where / is defined in the following. 
(1) Integral construction with unflanged stiifeners, 
J = (1/3)kt,? (21) 
(2) Z-stiffener panels 
J = (13) (h+ dp + da)t,? (22) 
Substituting the value of J in Eq. (20), and using E = 
2(1 + »)G with v = 0.3, we have 


= 0.715 (b (t t,)? (23) 


“Compare with Eq. (14) of reference 14. 
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Fic. 5. Relation between buckling parameter and aspect ratio, 


db. 


where //] takes the values and + + in in- 
tegral construction and Z-stiffener panels, respectively. 

For any given configuration, u is computed and the 
corresponding value of the buckling parameter y is ob- 
tained from the chart in Fig. 5. As stated before the 
critical stress can be obtained by using Eq. (14a). The 
value of the buckling parameter may vary from y = 
6.28 to y = 8.3 corresponding to the variation of u from 
infinity to zero. 

The stresses at these two extreme values of y corre- 
spond to buckling stresses of an infinitely long plate 
having its unloaded edges simply supported or fully 
fixed. 


MIXED INSTABILITY 


It is well known that many stiffened panels which 
suffer early plate buckling in short wave generally do 
not fail immediately, but are capable of carrying much 
higher loads. In panels of small slenderness ratio, the 
mode of buckling at failure may be Euler, in which the 
strong stiffeners with an effective width of skin act as 
columns. 

However, for practical design, we are interested in 
Z-stiffener panels having a large slenderness ratio 
(say L p 2 42). Many such panels are commonly 
employed in the low wing loading range [P L < 220 
(Ib. in. in), P = load per unit width]. The post- 
buckling behavior of such panels is governed by large 
deflections, and the nature of the problem is undoubt- 
edly complex. However, an approximate solution to 
this problem is suggested here. 

In test data’ '© many instances of panels with early 
skin buckling are found in which the average strain at 
failure is either coincident or nearly equal to 6, £ (é 
average failing stress). This suggests that final failure 
of panels depends on the failure of stiffeners. There- 
fore, in the panels, it is assumed that the skin, having 
suffered short wave buckling though undergoing in- 
creased deflection, is still capable of supporting the 
stiffeners, as far as the long wave buckling is concerned. 
Hence it is considered that a short wave buckled skin, 
and stiffeners, which have suffered negligibly small rota- 
tion at the instability, behave as a composite structure 
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and fail ultimately in long wave fashion. It is further 
assumed that the average failing stress in the skin is 
nearly the same as the average stress in the stiffeners. 
Therefore, the average failing stress is assumed to be 
the same as the critical stress at long wave buckling. 

On account of the two buckling modes being present 
in the same panel, the buckling phenomenon is referred 
to here as mixed (or double) instability. 


COMPARISON OF THEORY WITH EXPERIMENT 


Comparison of the theory is made with experiments 
covering panels reinforced by Z-stiffeners in the range 
0.25 t,/t 1.0, 15 b/t 75, and 12 h/t, 40; 
and integral construction with unflanged stiffeners in 
the range 1.4 < 2.0, 25 < b/t 35, and 9 < 
h/t, < 14. All the three types of instability that have 
been referred to here are examined. In Table 1 the 
result of analysis with respect to long wave buckling is 
given. In tests done by the NACA, the panels had 
flat ends, and had a fixity coefficient, ¢ = 3.75. How- 
ever, for computing the critical stress, only one wave 
is assumed between the two flat ends of the panel. In 
the tests done at Bristol Aircraft Ltd., on integral con- 
struction, the ends were clamped. 

Since long wave buckling, on account of the nature of 
the deformation and wavelength, is expected to cause 
immediate failure, the theoretical critical stress is com- 
pared with the average failing stress obtained in the 
tests. This mode of buckling is very important with 
respect to practical design, and, therefore, it is hoped 
that results given in Table 1, covering very wide range 
of configurations, will serve as a useful reference and a 
guide in design. 

Table 2 consists of results on short wave skin buck- 
ling. The theoretical critical stress is compared with 
that of the experimental critical stress, which was ob- 
tained as the stress at the load at which a plot of strain 
near the crest of the buckle first showed a decreasing 
strain with increasing load. 

In Table 3 the results of analysis on mixed instability 
are given, in which short wave buckling occurs as an 
initial instability and long wave buckling occurs at 
failure. The critical stress corresponding to long wave 
buckling is, therefore, compared with the average 
failing stress. In this table, the shortening per unit 
length at failure of the panel is also given. This was 
obtained in tests by the NACA as the average of strains 
indicated by four 6!» in. electrical resistance wire strain 
gages mounted on the quarter points along the length 
of two different stiffeners. 

Comparison of this with theoretically obtained value 
of o,,/E shows that the approximate solution suggested 
here for postbuckling behavior of the panels given in 
Table 3 is justified. 


CONCLUSIONS 


It is hoped that this article has further clarified the 
buckling phenomena of stiffened panels. In particu- 
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lar, the present analysis shows that buckling in the two 
wavelengths as given here is important in the design of 
compression panels. 

The two primary buckling modes, even though re- 
sembling Euler and local instability in wavelength, are 
distinctly different from them both in the type of defor. 
mation and in the critical stress. 

An extensive comparison of the theory with tests 
covering a wide range of both structural configuration 
and wing loading shows good agreement. 
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Structural R Ch teristi f 
ign of ome ructura esponse aracteristicsS Of a 
e e e 
hr arge Flexible Swept-Wing Airplane in 
1, are 
Rough Air 
tests 
THOMAS L. COLEMAN,* HAROLD N. MURROW,* ann HARRY PRESS** 
. 
Langley Aeronautical Laboratory, NACA 
mar) SUMMARY obtain a frequency decomposition of the strain responses 
a The paper presents an analysis of test measurements of airplane and thus to permit the determination of the contribu 
Thi motion and strain responses for a large swept-wing airplane in tions of the various airplane rigid-body and elastic- 
males flight through rough air. Power spectral methods are employed body motions to the strain responses at the various 
1949. in order to provide an indication of the frequency composition of wing stations. In addition, the various airplane rigid- 
the various responses and, in this way, to determine the contribu- m 2 ‘ a 
nsta- : é ays : ‘ body and elastic-body motions are examined in detail 
th tion of the various airplane rigid-body and elastic-body modes to 
' . the strain responses. The airplane rigid-body motions (pitch and in order to establish their contributions to the aero- 
7 vertical velocity) and elastic-body motions (deflections and dynamic forces. The results obtained in defining the 
sil twists) are examined in detail in order to assess the magnitude of contributions of the various airplane rigid-body and 
the aerodynamic forces arising from these motions. flexible-body modes to the strain responses and to the 
Pro. aerodynamic forces serve to provide a guide to the de- 
gree of refinement required in strain calculations. 
mit b = airplane span, ft. 
Vol. ( = mean aerodynamic chord AIRPLANE, INSTRUMENTATION. AND TESTS 
f = frequency, cps 
=m g = acceleration due to gravity Fig. 1 shows a plan-form view of the test airplane and 
V = airplane speed, the locations of some of the pertinent instrumentation. 
o(f) = power spectral density (see reference 5) 
stol lhe small x signs denote the locations of electrical wire 
= phase angle by which response lags gust input, deg. 
ep. ¢jne = root mean square value of flexible airplane response resistance strain gages installed on the front and rear 
o,; = root mean square value for quasi-static airplane refer- spars at four spanwise stations on the left wing and at 
ger ence condition the root station on the right wing. The approximate 
@ = increment in local or wing angle of attack, rad. spanwise location of the gages is given in terms of the 
= absolute value of complex number 
la- relative semispan station y (b 2) measured perpendic- 
of Subscripts used with a designate angle of attack change due to: ular to the fuselage plane of symmetry. The strain- 
29, B = wing bending velocity gage outputs are used here only as local strain indica- 
g = vertical gust velocity tions, and slow pull-up maneuvers are used for calibra- 
w- w = airplane vertical velocity 
a : tion purposes. The filled-in circles show the locations 
1 = wing streamwise twist i 
m | @ = pitch angle of optigraph targets installed near the front and rear 
1 r = rigid airplane spars at three spanwise stations. Deflections of the 
wing at these six points were recorded within the opti- 
wi INTRODUCTION graph mount on top of the fuselage. An NACA air- 
damped accelerometer (shown by the filled-in square 
Is p STRUCTURAL RESPONSE of airplanes in continu- 
ui ous rough air has been, and is, of continuing con- 
cern. In order to investigate this problem for contem- x STRAIN GAGES ff! 
5 porary airplanes, a flight investigation utilizing the © OPTIGRAPH TARGETS TF iy 
swept-wing B-47A airplane was undertaken by the fe 
- NACA in 1955. Some of the early results from this 
inves > 2e ‘ 2 tere 
investigation hav e been published in references | to 
rhe purpose of this paper is to present a more detailed VANE 
analysis of the test data in order to describe the signifi- 
; cant factors affecting the strain responses, and, particu- 
larly, the effects of airplane dynamic flexibility. Power 
spectral methods of analysis are applied in order to. 
Presented at the Aeroelasticity—I Session, Twenty-Sixth 
Annual Meeting, IAS, New York, January 27-30, 1958. 
* Aeronautical Research Engineer. 
** Supervisory Aeronautical Research Engineer. Fic. 1. Test airplane and instrumentation 
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symbol) was installed near the center of gravity to 
measure normal acceleration. In addition to the 
accelerometer shown at the center of gravity, acceler- 
ometers were installed at 22 other locations to measure 
local accelerations at various locations on the wing and 
fuselage. A mass-balanced metal flo v-direction vane 
was mounted on the nose boom in order to obtain a 
measure of the angle of attack of the airplane for use in 
determining the vertical gust velocities. Measurements 
also were made of the control movements and of the 
angular displacements and velocities at the center of 
gravity about the three axes. 

The basic test measurements to be used were ob- 
tained during test runs of several minutes duration in 
clear-air turbulence at pressure altitudes of approxi- 
mately 5,000 ft. and 35,000 ft. Most of the results to 
be presented are for the low-altitude test condition. 
The flight-test conditions and the duration of the data 
sample for the two altitudes are summarized in Table 1. 
It should be mentioned that, for this airplane, most of 
the airplane mass is concentrated in the fuselage and, 
for the present investigation, the ratio of the wing 
weight to total weight was approximately 0.36. The 
flights in rough air were made by using a minimum of 
control motion—that is, minor deviations of the air- 
plane from the prescribed altitude and heading were not 
corrected by the pilot, and major deviations were cor- 
rected only by gradual control movements. The 
flights were made without the use of a yaw damper. 
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g taken during these pull-up maneuvers are used to 
determine quasi-static reference strains for comparison 
with the strain measurements taken during rough-air 
flight, as will be described. 


CHARACTERISTICS OF STRAIN RESPONSES 


Time-History Measurements 


In order to illustrate the general characteristics of the 
various airplane responses to the rough air, sample time 
histories of some of the measurements are shown jn 
Fig. 2. The upper curve represents the time history of 
the vertical gust velocity as determined from measure- 
ments of the angle of attack vane, the pitch, and the 


airplane vertical velocity as described in reference 4, | 


The next two time histories are bending strain responses 
as measured at an inboard station and a midspan sta- 
tion. Also shown in the Figure are the time histories 
for the c.g. normal acceleration, airplane vertical 
velocity, pitch angle, roll velocity, and wing deflection 
at an outboard station. 

Examination of the various time-history measure- 
ments indicates that both the gust velocity and the 
various airplane responses are erratic and irregular in 
character. This characteristic of the rough air and the 
airplane responses to the rough air is general and is the 
basic reason for the recent applications of random- 
process techniques in the analysis of these problems. A 
second point of interest is the wide differences in ire- 
quency content between the various responses. The 
angular motion responses of both pitch and roll velocity 
are predominantly long wave in character, while the 
gust velocity, normal acceleration, and strain responses 
reflect more of the higher frequencies. 

Considering the time histories in Fig. 2 in further de- 
tail, a comparison of the magnitudes of the angle of 
attack changes induced by some of the motions with the 
direct gust-induced angle of attack is of interest. Fora 
test airspeed of 700 fps, the fluctuations of +10 fps in 
the gust-velocity time history correspond to gust- 
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POWER SPECTRAL 


In addition to the flights in rough air, slow pull-up 
8 us DENSITY, &(f) 
maneuvers were also performed in smooth air at a 
variety of test conditions in order to determine the 
strains, deflections, and twist per unit load under quasi- 
static conditions. The measurements of the strain per 
TABLE 1 
e.g. Duration ol 
altitude, Weight, Mach per cent sample, ie) | 2 3 o S 6 
ft. Ib. Number min. f, 
5,000 113,000 0.63 20.0 4.0 
MIG. yer spectri -nsity -ar-Spd strains for 
35,000 112,000 0) 64 20 3 15 Fic. 3. Power spectral density of rear-spar bending strains fo 
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STRUCTURAL 


induced angle of attack changes of the order of +1°. 
The airplane angle of attack changes arising from both 
the airplane vertical velocity and the pitch angle can be 
readily seen to be of the same order of magnitude. The 
roll-velocity motions for this airplane (with a 57-ft. 
semispan) also introduce angle of attack changes at the 
wing tips of almost the same magnitude. 

Another point worth noting is that, whereas the 
general patterns of the various responses are similar, 
significant differences in the phases exist for some of the 
peaks, as may be seen by comparing the responses at 
time f = 3.5 sec., shown by the dashed line in the Figure. 
It may also be noted that, although the strain peaks and 
deflections are in phase, they both lag the gust and c.g. 
acceleration peaks. The pitch angle is likewise out of 
phase. These phase relations, as will be seen, play an 
important part in the airplane response. 


Power Spectral Analysis of Bending and Shear Strains 


Basic Method——One of the principal aims of the pres- 
ent analysis is to provide a measure of the effects of 
dynamic flexibility on the strains. For this purpose, a 
basic comparison of interest that may be obtained from 
the flight test data is a comparison of the measured 
strains with the strains that would be developed if the 
same aerodynamic loads were applied in a static man- 
ner. The particular method of arriving at this 
dynamic-static comparison employed in the present 
paper consists of the following steps. 


(a) The strain measurements in rough air are con- 
verted to equivalent quasi-static system accelerations 
by dividing the turbulence induced strains at various 
stations by the strain per g measured at the correspond- 
ing station in slow pull-up maneuvers. The results of 
the analyses are spectra of quasi-steady system accelera- 
tions required to give the spectra of the strains measured 
in turbulence. 

(b) An approximate power spectrum of the airplane 
system acceleration was obtained by analyzing airplane 
¢.g. accelerations (results reported in reference 1 sho v 
that this is a reasonable approximation). It can be seen 
that the comparison of this power spectrum with those 
obtained by the procedure in (a) above will provide an 
indication of the structural dynamic effects. 

The differences between the loads obtained from the 
strains and the measured loadings provide a direct 
measure of the dynamic effects on the strain response. 
This procedure has the particular advantage of yielding 
a common quasi-static reference loading condition for 
all the strain stations. 

Bending Strain —The power spectra of the bending 
strain on the rear spar at the four spanwise locations are 
given in Fig. 3. The common quasi-static reference 
power spectrum, which, as indicated earlier, was based 
on the c.g. accelerations, is also shown in the Figure. 
The difference between this reference curve and the 
Spectra at the various stations provides the desired 
measure of the effects of dynamic flexibility on the 
Strains at the various frequencies. Inasmuch as these 
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Power spectral density of rear-spar shear strains for four 
spanwise stations. 
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power spectra are, in effect, spectra of equivalent 
accelerations, the units for the ordinate scale are 
(ft. sec.*)?/eps. 

Examination of Fig. 3 indicates that almost all of the 
power at the various stations is concentrated at fre- 
quencies below 2 cps. In this frequency region, the 
pover appears to be concentrated in three lobes or 
peaks —one at about 0 to 0.2 eps, one at 0.5 eps, and the 
last at 1.5 eps. The peak in power at 0 to 0.2 eps can be 
seen to increase rapidly for the more outboard positions 
and is a reflection of the rolling motions associated with 
the airplane Dutch roll mode and also, to some extent, 
the result of pilot aileron control motions. Thus, this 
power peak has no direct bearing on dynamic aero- 
elastic effects. The peak in power at 0.5 eps is of about 
the same magnitude at all stations and is a reflection of 
the strains associated with the airplane short-period 
mode. The peak in power at 1.5 epsis a reflection of the 
fundamental wing-bending mode. This peak is seen to 
be smallest for the root station and to increase progres- 
sively toward the midspan stations. In addition, at 
higher frequencies, evidence of the effects of higher 
structural modes is apparent at the outboard stations. 

Comparison of the various strain spectra with the 
quasi-static reference curve provides an indication of 
the dynamic flexibility effects. If the differences at 
very low frequencies (which, as indicated, are associated 
with rolling motions) are neglected, flexibility can be 
seen to introduce large dynamic strain amplifications at 
frequencies above 1 cps, particularly in the neighbor- 
hood of the fundamental bending mode at 1.5 eps. At 
higher frequencies, the amplifications are large but the 
total power is relatively small. 

In order to show the overall effect of the airplane 
flexibility on the strains, the ratios of the root mean 
square strains to the root mean square strain for the 
reference spectrum (¢y;,: ¢,,;) for the various stations 
are also shown in Fig. 3. In computing these values, 
the incremental power at f < 0.2 eps associated with 


rolling motions was not used. The strain ratios vary 
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total bending and shear strains [y/(b/2) = 0.08, rear spar]. 


from 1.12 to 1.66 and indicate that the strain amplifica- 
tions are smallest at the root station and are largest at 
the 0.60 semispan station. 

It is also of interest to consider separately the effects 
on the root mean square strains of the contributions due 
to the lateral motions and the higher structural modes. 
From the power spectra in Fig. 3, it can be seen that for 
the root station, the effects of the lateral motions and 
the higher structural modes on the root mean square 
strains appear to be negligible. For the outboard sta- 
tion, the effects of both the lateral motions and the 
higher structural modes appear to be significant, and 
each acts so as to increase the root mean square values 
by roughly 5 to 10 per cent. 

A similar analysis was made of the front-spar bending 
strains. Although these results are not given, the 
strains had essentially the same spanwise variations as 
were found for the rear spar. 

Shear Strains—The power spectra of the shear strains 
at the four spanwise wing positions are shown in Fig. 4 
together with the quasi-static reference power spectrum. 
The shear-strain power spectra exhibit the same general 
characteristics as were noted in Fig. 3 for the bending 
strains. However, two differences between the shear 
and bending strain responses should be noted. First, at 
low frequencies the effects of the Dutch roll and aileron 
control motions appear greater on the shear strains than 
was noted on the bending strains. Again, the effect of 
the lateral motions progressively increases from the 
root station to the outboard station. Secondly, the 
effects of the higher elastic modes are more pronounced 
in the case of the shear strains than in the case of the 
bending strains. Examination of the curves at fre- 
quencies above 2 eps indicates the presence of several 
elastic modes that were not noticeable in the bending 
strain power spectra. For example, strain peaks associ- 
ated with elastic modes at approximately 2.2 cps, 5.2 
cps, +.5 eps, and 5 eps are evident, although the effects 
of these modes vary between wing stations. 

The dynamic amplifications of the shear strains 
appear to result primarily from the first wing-bending 
mode, except at the root station, where little evidence 
of this mode is noticed. The overall effects of wing 


flexibility on the shear strains are also indicated by the 
ratios of root mean square strain for the flexible air- 
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plane to the root mean square strain for the quasi-static 
reference airplane. Again, the incremental power at 
low frequencies associated with the rolling motions jg 
omitted in determining the strain ratios. The strain 
ratios range from 0.94 to 1.72 and are similar in magni. 
tude to the ratios shown for the bending strains. 

It is of interest to consider separately the effects on 


the root mean square shear strains of the contributions 


of the lateral motions and the higher structural modes, 


For the root stations, the contributions of both of these | 


sources to the root mean square values appear negligible. 
However, for the outboard stations, the lateral motions 
and the higher structural modes increase the root mean 
square strains by approximately 20 per cent and 10 per 
cent, respectively. 

The effects of flexibility on the front-spar strains 
were of the same order as that shown for the rear spar. 
However, the front-spar strains showed a less orderly 
spanwise variation than the rear-spar strains, suggesting 
the need for additional study of torque strains in rough 
air. 

Effect of Antisymmetric \lodes—The power spectra of 
bending and shear strains presented in Figs. 3 and 4 
reflect the strains arising from both symmetrical and 
antisymmetrical airplane and elastic response modes. 
In order to examine the contribution of the antisym- 
metric modes to the strains, a time history of the dif- 
ference between the strain time histories at the root 
station on the leit and right wings was obtained. The 
power spectra of the antisymmetric contributions to 
the bending and shear strains are given in Fig. 5 for the 
rear spar at the root station. (This spectrum is the 
spectrum of one-half the difference in strain between 
the right and left wing.) For comparison, the spectra 
of the total bending and shear strains measured at this 
station are also given. 

The spectrum of the bending strain due to the anti- 
symmetric modes shows power peaks at several fre- 
quencies. The strain peaks at very low frequency (be- 
low 0.25 cps) show the strain power associated with 
rolling motions and aileron control movements. An- 
other power peak exists at a frequency of approximately 
2.2 eps and is a reflection of an antisymmetric elastic 
mode. The shear-strain power spectra shows the same 
general pattern but, in addition, indicates the presence 
of a significant antisymmetric mode at a frequency o/ 
4.5 eps. 

In comparison with the spectra of the total bending 
and shear strain at the root station, the strain assoct- 
ated with the antisymmetric modes is negligible in the 
frequency region below 1.5 cps. At higher frequencies, 
however, the antisymmetric modes appear to contribute 
significantly to the total strain response. Unfortu- 
nately, similar comparisons could not be made for the 

outboard stations. 

Variation With Altitude—tIn order to illustrate the 
variation in the characteristics of the strain responses 
with altitude, the power spectra of bending strain meas- 
ured at two altitudes (5,000 ft. and 35,000 ft.) are given 
in Fig. 6 for the rear spar at the root station. The power 
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spectra of the quasi-static airplane reference condition 
are also shown in the Figure. As previously mentioned, 
the airplane weight and Mach Number were approxi- 
mately the same for the tests at the two altitudes. 

The difference between the levels of the strain spectra 
in Fig. 6 is primarily a reflection of lower dynamic 
pressure for the high-altitude tests, since the intensities 
of the turbulence at the two altitudes were approxi- 
mately the same in terms of true gust velocities. Com- 
parison of the strain spectra with the reference spectra 
in Fig. 6 shows that the strain amplification associated 
with the fundamental bending modes is relatively much 
greater at 35,000 ft. than at 5,000 ft. This increase in 
the excitation of the elastic modes is a reflection of the 
decrease in the aerodynamic damping at the higher 
altitude. The overall strain amplification is shown by 
the ratios of the root mean square strains Ofex ref 
noted in the Figure for the two altitudes. Comparison 
of these two ratios shows that the strain amplification 1s 
approximately 15 per cent greater at the high altitude. 
Similar results were also obtained for the other strain 
stations. 

Implications—The foregoing analysis of the strain 
power spectra appears to have a number of implications 
in regard to structural design calculations. First, it is 
clear that for airplanes of the present type, the effects of 
dynamic flexibility in rough air significantly increase 
the wing strains over most of the airplane span. The 
dynamic amplifications vary widely with spanwise 
position, and, thus, the distribution of strain across the 
airplane span in rough air is quite different for the gust 
condition than for the maneuver condition. As a con- 
sequence, comparisons of gust and maneuver loading 
conditions have to be considered in greater detail, with 
particular consideration given to the consequences of 
the dynamic-gust-response effects. In regard to the 
degree of refinement required in gust-response analyses, 
it would appear that for the root stations, analyses con- 
sidering the degrees of freedom of airplane vertical 
motion, pitch, and the first wing bending mode might 
be expected to yield reliable estimates of the strain. 
For the outboard stations, however, the contribution to 
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Fic. 6. Power spectral density of rear-spar bending strain for 
two altitudes [y/(b/2) = 0.08]. 
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Fic. 7. Frequency-response functions of angle of attack change 


due to airplane vertical velocity and pitch for vertical gusts. 


the strain arising from the airplane lateral dynamics and 
the higher flexible modes, both symmetric and antisym- 
metric, may also have to be considered, depending upon 
the accuracy desired. 


CONTRIBUTIONS OF AIRPLANE MOTIONS TO 
AERODYNAMIC FORCES 


In the preceding analysis, it was indicated that, at 
low frequency, the strain responses are principally asso- 
ciated with the airplane longitudinal short-period mode. 
At the higher frequencies the airplane vibration modes, 
particularly the first wing-bending mode, contribute 
substantial amplifications to the strain responses. The 
strain responses in both these frequency regions may 
be expected to depend in large part upon the magnitude 
of the aerodynamic forces which arise from the rigid- 
body airplane motions in the short-period frequency 
region and from the elastic-body motions in the first- 
bending-mode frequency region. In the remainder of 
this paper, an effort is made to assess the contributions 
of the rigid-body and elastic-body motions to the aero- 
dynamic forces. The basic approach to be used is to 
examine the manner in which the various motions 
affect the wing or local angle of attack. 

Rigid-Body Motions—In order to assess the effect of 
the airplane pitch and vertical velocity on the aero- 
dynamic forces, frequency-response functions of the 
angle of attack due to pitch @ and vertical velocity 
(w WV’) to continuous rough air were determined from 
the test data. The frequency-response functions were 
obtained by the cross-spectrum method described in 
reference 5, and the results obtained are shown in Fig. 
7. These curves describe the amplitude and phase of 
the response to unit sinusoidal gust velocities at the 
various frequencies. For present purposes, the fre- 
quency-response functions are given in terms of the 
angle of attack increments induced by the two responses 
(nose-up pitch and downward airplane vertical velocity 
are taken as positive). In determining these frequency- 
response functions, the vertical gust velocities derived 
from the vane angle of attack measurements were used 
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Fic. 9. Frequency-response functions of angle of attack 
change due to response of streamwise twist and wing-bending 
velocity to vertical gusts [y/(b/2) = 0.74]. 
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as the gust input, and the airplane vertical velocity was 
determined from the normal acceleration measurements 
at the center of gravity. 

Inspection of Fig. 7 shows that the amplitudes of the 
angle of attack have a similar variation with frequency 
and that the maximum responses occur at about 0.5 
cps. In comparison with the angle of attack increment 
due to a unit gust velocity (indicated in the Figure by 
the arrow on the ordinate), the magnitudes of the pitch 
and vertical-velocity responses appear to be significant 
up to frequencies of about 2 cps. 

The effect of the pitch and vertical-velocity responses 
on the net airplane angle of attack change depends on 
the phasing of the responses relative to the gust veloc- 
ity. Consideration of the phase angles in Fig. 7 shows 
that at low frequency the angle of attack arising from 
both responses is nearly 180° out of phase with the gust, 
and thus an alleviation of the gust-induced angle of 
attack is realized. However, as the frequency increases, 
the phase relations for both motion responses change 
and, in fact, become in phase with the gust at some fre- 
quencies. Consequently, at these frequencies, the 
motions will tend to amplify the gust-induced angle oi 
attack. The net effects of the motion responses on the 
angle of attack are summarized in Fig. 8. 

The ordinate in Fig. 8 represents the ratio of the 
amplitudes of the net angle of attack change to the 
direct gust-induced angle of attack change. The three 
curves represent the effects of the addition to the gust 
angle of attack a, of the angle of attack due to: (a) 
pitch motions, a,; (b) airplane vertical velocity mo- 
tions, a,; and (ec) both pitch and vertical motions, 
ay + ay. These curves were obtained by vectorial addi- 
tion of the frequency-response-function vectors, as 
indicated by the sketch on the Figure, for f = 0.50 eps. 
They represent the ratios of the amplitudes of the air- 
plane angle of attack change (with the effect of the 
specified motion added) to the direct gust-induced 
angle of attack. These curves can thus be considered 
to represent an angle of attack amplification function 
at the various frequencies. 

From consideration of the results shown in Fig. 5, it 
is apparent that, at low frequencies, both airplane verti- 
cal velocity and the pitching motions alleviate the 
angle of attack induced by the gusts. However, at 
frequencies above about 0.5 cps to almost 1.5 eps, the 

situation changes drastically with both the pitch and 
vertical velocity motions tending to amplify the air- 
plane angle of attack. At higher frequencies, both 
motions are negligible. Because of the phase relations, 
the combined effect of the pitch and vertical velocity is 
much the same as the result for pitch alone, reducing 
the angle of attack at low frequencies but amplifying 
the angle of attack between frequencies of 0.3 to 1.25 


cps. 

In order to determine the overall effects of these rigid- 
body motions on the overall angle of attack and loads 
in rough air, the frequency spectrum of the gust input 
must also be considered. It has been estimated on this 
basis that the overall airplane loads are decreased by 
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about 10 per cent as a result of airplane vertical motions 
but are increased by about 35 per cent as a result of the 
pitching motions. 

Elastic-Body Motions—In order to examine the 
effects of the elastic-body motions on the local angle of 
attack, the frequency-response functions for the angle 
of attack due to the wing streamwise twist and the 
wing vertical velocity due to bending alone were also 
determined from the optigraph records and are shown 
in Fig. 9. These results are for the 74 per cent semi- 
span station. Again the frequency response functions 
are given in terms of the angle of attack change per unit 
gus velocity. (For this purpose, twist of the wing 
leading edge upward yields a positive angle of attack 
change, and wing bending velocity downward yields a 
positive angle of attack change.) The magnitude of the 
angle of attack changes due to bending and twist can 
be seen to be small at the low frequencies relative to the 
direct gust-induced angle of attack (indicated by the 
arrow on the ordinate). As the frequency increases, the 
angle of attack changes due to both flexible motions 
become very large, particularly in the neighborhood of 
the first bending mode, and are, in fact, at these fre- 
quencies several times larger than the gust-induced 
angle of attack or the angle of attack changes due to 
rigid-body motions as shown by Fig. 7. The phase 
behavior shown by the lower curves is also significant. 
As indicated on the power part of the Figure, at low 
frequencies, the twist lags the bending velocity angle by 
90°. At the higher frequencies, the phase relations 
change and the lag is increased to about 180°. 

Fig. 10 shows the net effects of these flexible motions 
on the local angle of attack. The Figure gives the ratio 
of amplitudes for the net angle of attack (when the 
flexible motions are included) to the angle of attack for 
the rigid airplane (as obtained from Fig. 7). As before, 
the ordinate can thus be considered an angle of attack 
amplification ratio and in this case represents the 
amplification due to flexible motions. Three curves are 
shown to represent the effects on the local angle of 
attack of: (a) wing-bending velocity alone, ag; (b) 
twist, a7; and (c) both bending velocity and twist, 
aj + ar. These curves were also obtained by vectorial 
addition, as indicated by the sketch in the upper part 
of the Figure. 

Examination of the results in Fig. 10 indicates that, 
at low frequencies, only the twist is important and tends 
to alleviate the gust-induced angle of attack, as would 
be expected for a swept-wing airplane. At the higher 
frequencies, both motions increase the angle of attack 
by considerable amounts. However, the 180° out-of- 
phase character of the bending velocity and twist 
causes the angle of attack changes due to both motions 
to be reduced below that due to twist alone, as indicated 
by the sketch for f = 1.5 eps. The overall angle of 
attack changes at these frequencies are still large and 
are close to three times those obtained for the rigid 
body. (At these higher frequencies, the rigid-body 


motions are small, as indicated by Fig. 7, and the rigid- 
body angle of attack change is very close to the direct 
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gust-induced angle of attack change.) The effects of 
these large angle of attack changes are, of course, some- 
what moderated by the attenuation of the magnitude 
of the lift for oscillatory motions at high frequencies. 
Nevertheless, the amplifications of the aerodynamic 
loading can still be expected to be large. 

Implications—The preceding results on the contri- 
butions of the longitudinal and elastic motions to the 
aerodynamic forces in rough air appear to have two 
major implications in regard to strain calculations. 
First, it appears that a good representation of the longi- 
tudinal motions, particularly the pitching motions, is 
required in order to define adequately the aerodynamic 
forces in the lower frequency range (below approxi- 
mately 1 cps). Second, it appears important to define 
the elastic motions (wing bending and twist) accurately, 
inasmuch as these motions are the dominant sources of 
the aerodynamic forces in the frequency region of the 
first wing-bending mode. 

CONCLUSIONS 

An analysis of flight test data obtained with a large 
swept-wing airplane in continuous rough air has indi- 
cated the following results. 

Significant contributions to the wing strain responses, 
both bending and shear at various spanwise locations, 
arise from rolling motions associated with the lateral 
Dutch roll mode, longitudinal short-period motions, and 
elastic modes of vibration. 

The effects of airplane flexibility on the strain re- 
sponses are predominantly associated with the first 
wing-bending mode. At the outboard stations, signifi- 
cant contributions also arise from higher frequency 
modes, both of a symmetric and antisymmetric nature. 

The effects of airplane flexibility cause a moderate 
amplification in the wing strains at the root (10 to 20 
per cent) but a much larger strain amplification (50 to 
70 per cent) for the midspan region. 

For flights at the same Mach Number, the effects of 
airplane flexibility cause substantially greater strain 
amplifications at high altitudes than at low altitudes, 
presumably resulting from the lower dynamic pressure 
and the reduced aerodynamic damping. 

The effects of the airplane rigid-body motions, and 
particularly the pitching motions, yield large amplifica- 
tions in the wing angle of attack, especially at frequen- 
cies in the neighborhood of the short-period mode. 
These amplifications of angle of attack result in signifi- 
cant increases in the aerodynamic loading at these 
frequencies. 

The effects of the airplane elastic motions, both wing- 
bending velocity and streamwise twist, have significant 
effects on local angle of attack variations on the wing 
and the associated aerodynamic loading. At low fre- 
quencies (< 0.8 cps), favorable streamwise twists asso- 
ciated with the swept wing tend to reduce local angles of 
attack and alleviate the aerodynamic loading. At the 
higher frequencies (1.5 cps) these vibratory motions, 
however, tend to increase the angle of attack changes 
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A Means of Improving the Static Performance 
of Cruise-Designed Shrouded Propellers | 


Aeronautical Research Engineer, Aerodynamics Laboratory, 
David Taylor Model Basin, Washington, D.C. 


March 28, 1958 


O~ OF THE FEATURES of shrouded propellers that make them | 
attractive for use in certain types of VTOL aircraft is their 
capability of producing high static thrust. However, a problem 
in applying them to this use arises from the requirement of VTOL 
aircraft to operate effectively in two regimes—static flight and 
cruising flight. Shrouds designed for the static condition have 
high drag at cruising speeds, while cruise-designed shrouds in 
static operation experience losses, due largely to separation, that 
reduce shroud thrust. Therefore, it is highly desirable to devise 
shrouds that can be easily converted from one configuration to 
another in flight. One such scheme is a shroud designed for 
cruising conditions but equipped with a leading-edge slat that 
can be extended to improve its static performance (Fig. 1 
Static tests of this shroud have been conducted in order to deter- 
mine the improvement afforded by the slat. 

Shroud effectiveness is defined as the ratio of the shroud thrust 
actually developed to the maximum value theoretically possible 
One-dimensional momentum theory states that for nondiffusing 
shrouds under static conditions shroud thrust 7, is equal to pro- 
peller thrust 7, when there are no shroud losses (see, for ex- 
ample, reference 1). Consequently, shroud effectiveness may be 
expressed as 7,/T,. The value of shroud effectiveness approaches 
unity as shroud losses are reduced. 

In the experimental investigation herein reported, 
thrust and propeller thrust were measured with various slat ex- 
tensions. The variation of shroud effectiveness with ratio of slat 
extension / to shroud diameter D is shown in Fig. 2. At full 
practicable extension (//D 


shroud 


= 0.325), the improvement in 7,/T, | 
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Fic. 1. Slatted-shroud model. 
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Fic. 2. Effect of slat extension ratio on shroud effectiveness. 


was about 100 per cent as compared with the cruise configuration 
= 0). Note that for //D = ©—i.e., with the slat entirely 
removed-—ideal shroud effectiveness is very nearly realized. 

It is expected that similar increases in shroud thrust would be 
realized at much smaller slat extensions if the slat shape were 
modified to produce larger openings at given extensions. 
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On One-Dimensional Temperature 
Distribution in Two-Layered Slabs With 
Contact Resistance at the Plane of Contact 


Paul Seide 


Space Technology Laboratories, a Division of The 
Ramo- Wooldridge Corporation, Los Angeles, Calif. 


April 1, 1958 

sang YEARS AGO in reference 1 the writer presented a solu- 
tion of the problem of the one-dimensional transient tem- 

perature distribution in a two-layered slab (or two insulated rods) 
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with contact resistance between the two layers when one surface 
is subjected to a sudden rise in temperature, 7), and the other 
end is held at the initial zero temperature of the slab 
therein, the solution is limited to the case where the thicknesses, 
h,, and diffusivities, a, of the two layers satisfy the relationship 


As derived 


h,? a; = he? 


The present note gives the general solution of the problem and 
simplifies the result of reference 1. 

We denote the conductivities of the layers by A, and Ko», the 
temperatures in the layers by 7) and 7» (functions of distance x 


and time ¢), and the contact conductivity by / Then the 
temperatures are governed by the equations 
a, = O7;/0t, OL Sh, (2a 
= O72/01, Sh (2b 
subject to the conditions 
T, = T, atx, = 0, t>0 (3a 
= t>0 (3b) 
K,(07;/0x;) = = 
— atx; = = 0 (3c 
T, = T; = O for allx,¢ = 0 (3d) 


Consider first the steady-state condition in the slab for infinite 


time. For thiscase 7) and are given by 
(he/Ke) + (4a) 
= To [1 — + 
+ (1/U)|{ (4b) 
which satisfies equations (2a), (2b), (3a), (3b), and (3c For 
the transient state let us denote the solution by 
T; = T;,, + > e— Ant (x;) i= 1,2 (5) 
n 
where 
Xin(X1) = Ay, sin V (x /hy 
> 


Xon(X2) = Ao, sin V [1 — (x2 hs)|$ 
which satisfies equations (2a), (2b), (3a), and (3b To satisfy 


equation (3c) we must have 
Ay,(Ki/hy) V cos + 
A», (Ke/h2) V Brhe?/az cos V = (Fa) 
(11/2) (Ki/hy) V Cos V ay + 
Usin — (K2/he) V 
cos V + U sin V = (7b) 


We assume that h,?/a; ¥ he?/az. Then, from the condition of 


nontrivial A;, and A»,, the characteristic equation for 8,, is 


{tan V V 


+ (1/U) =0 (8a) 
i=1 
and a relationship between the coefficients is 
hy) Vv B, hy? COS a; = 
— Ao,(Ke/he) V cos V By a2 = An (Sb) 


From equations (3d), (5), (6), and (8b), the conditions for the 


determination of A,, are 


A, {sin hy? ay /(K1/hy) V cos 


V + = 0 (9a) 
> An V By [1 — (x2 hy)| /(Ke/he) X 
n 

V B,he?/a2 cos V — = 0 (9b) 


Let us now multiply equation (9a) by 


(hy2/ay) [sin V V cos a] 


~ 
| 
| 
‘ ‘ 
‘oud | 
slat 
full 
| 
2 
| — 
| 


integrate over the region 0 < x,/h; < 1, and add the result to 
equation (9b) multiplied by 


2 2/ 1 
= B, hj? / aj cos V By hi?/ a: 
It ean be shown by direct integration and use of the characteristic 
equation (Sa) that the quantity in brackets vanishes if m + n, 
and is equal to 


/ 
En {sin a; 2vV B,h;? a;] 
where ¢, = (lla) 
Bn i=1 (K;/h,;) (A + cos 2V a;) 


if m = n, from which follows the relation 
An = —(1/¢n)Ty see V an (11b) 
If hy?/ = he? = y, the solution must be obtained somewhat 
differently. We note that the quantity cos V 8,y is a factor of 
equation (7a). Then we may have 
cos V By = 0 (12a) 


in which case, from equation (7b), 


Ain = Ao, = A, (12b) 
or we may have 
(Ay hy) Ain = — (Ks hy) Aon = B,, (12¢) 


in which case, from equation (7b), 
(tan V BaW/V Baw) + + (2/K2)|} = 0 (12d) 


The transient temperatures are given by 


7, + A,e~ sin [((2n — 1)/2] X 


n=1 


mixi/Iy) + 8" sin V (18a) 
n 
T= Ts, + A,e 2n+D/2] 
sin — 1)/2|r[1 — (x2/he)| 


8"! sin V Buy [1 — (x2/he)] (8b) 


n 


where 3, is defined by equation (12d). Since both (x,//,) and 
{1 — (x2/he)| vary from zero to unity, let us denote them, as in 
equation (10), by the quantity & The conditions for the deter- 
mination of A, and B, are then given by 


A, sin [(2n — + |B, /(Ki/h)| X 
n 


n=1 
sin — = —Th,, 
; (14) 
A, sin [(2n — 1)/2]ré — [B,/(K2/hz)| X 
n=1 n 
sin E = — Tr, 


from which we obtain 


A, sin [(2n — 1)/2]ré = 


‘m1)T + (K2/hs) + (Ky hy) (15a) 


DB, sin VB (To, — + (e/K2)| 


n 


We note that 


f 2n- 1 si 2n 1 jo ifmatn 16 
sin sin - = (lba) 
0 2 (1/2) ifm=n 


and, by use of equation (12d), 


V B,,hi?/a; cos V B,,h;?/ a: 
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(hy?/ az) [sin V (X2/h2)/V az cos V By hy? ay 


and integrated over the region 0 < 1 — (xe/ho) < 1.) We the 
obtain 


2 | 2 
sin V 8B,,h;2/a; de l= — T, 
Mm ay 


sin V sin VB,, EdE = 
0 


sO ifm +n 


(1/2) [1 — (sin V28,y/2V ifm =, 
Then, it is readily determined that 
An = + (K2/he)|} [47 + (17, 
Bn = —|T)/(2V Bu — sin 2V X 


(hy K;) (hy l7b 
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Approximate Formulas for Thermal-Stress 
Analysis 


D. J. Johns 
Lecturer, The College of Aeronautics, Cranfield, England 
April 3, 1958 


—— BASIS of any thermal-stress analysis is the determination i 
of the temperature distributions in the structure. For arbi- | 


trary flight histories, the determination of such distributions is 
rather tedious and not completely general. This latter fact 
handicaps optimization studies in the project design stage wher 
it is desirable to be able to express the thermal-stress distributions 
in a general manner. 

In this note, general expressions are derived for the thermal- 
stress distributions in a typical I-section using similar assump- 
tions to those of Biot! with the following additions and exceptions 


A = area 

E = Young's Modulus 

x,y = coordinates defined in Fig. 1 
a = coefficient of expansion 
stress at position y,x 

6y,2 = temperature at position y,x 


Suffixes F, W refer to flange and web, respectively. The flange 
and web materials are assumed to be dissimilar and their proper 
ties invariant with temperature. There is no temperaturé 
drop across the flange to web joint. 

There are two distinct phases in the heating of an I-section 


Initially, the temperature does not begin to rise at the midpoint | 


of the web and everything happens as if the web depth were 
infinite. 

During this first phase the web temperature distribution * 
approximated by 


6, =H {1 — 
6, = 0 for (J — y) > qi) 


where q; is a generalized coordinate and may be called the ‘pent 
tration depth.”’ The time ¢ = 4, at which g; = / (i-e., the tem 
perature at the web center begins to rise) is called the “transit 
time.”’ 

For the problem of a uniform slab with one face at constatt 
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temperature, 
= 3.36 V kut/c, (2) 
= 0.0885 (3) 
For the more realistic case of an I-section with 4 = f(t) the 


ibove results are not strictly correct but can be assumed with 
little error (reference 1, Fig. 5). 

During the second phase of heating (¢ > ¢,) the temperature 
't the center of the web is rising, and the web temperature dis- 


tribution is approximated by 
6y = O(y/l)? + — (y/l)?] (4) 
During all phases of heating, if the temperature of the flange 
remote from the web is 6, then the flange temperature distribu- 
tion is approximated by 
6, = [1 — (x/d)]? + @(a/d) [2 — (x/d)] (5) 
4, is easily determined for any arbitrary flight program since 
conduction of heat at element 2 may be ignored. If radiation 
losses are included, 9 must satisfy the differential equation 


acpO, = — 62) — (6) 


It only remains, therefore, to express temperatures 6; and 65 as 
functions of 02. 
Hoff? has given the following expression for the first phase of 
heating: 
= (1 — p*) & ( 


where 1/p* = 1 + 2(a/2a;) (ke/kw) (8) 


and this approximation may be assumed for ¢ < 4). 

Alternatively by extending Biot’s analysis of an I-beam using 
dissimilar flange and web materials, the differential equation de- 
fining #; may be approximated by using a modified heat capacity 
(ac)’ = B’ac in Eq. (6) where 


= 1+ (a)/a) (kw/ke) = 1 + (aie,/acp) (9) 


Although Biot’s original analysis strictly only applied for ¢ < f; 
he justified its use for all values of fin one application. 

To determine 6, Biot’s slab analysis will be assumed for the 
case when & and 6; are arbitrary functions of time. The resultant 


equation obtained is 


6; + 4.576; = 6, — 1.07516, (10) 
When 6; is constant, Biot obtained 
6; + 4.671)65 = 6; (11) 


The difference in coefficient of 6; is due to a numerical error in 
reference | 

Therefore, summarizing, 6 is found from an equation similar 
to Eq. (6). 
found in the literature. 


Numerous examples of such calculations may be 
To determine 6, a modified value of the 
heat capacity (ac)’ = B’ac probably gives the best solution, al- 
though, in the first heating phase Hoff’s results, Eq. (7), may be 
In the 
general case a numerical solution of Eq. (10) will be most con- 


more convenient. 6; may be approximated by Eq. (10). 
venient 

The general expression for the thermal stresses in a symme- 
trically heated I-section is 


Therefore, using the temperature distributions of Eqs. (1), (4), 
and (5) the following expression for oy,, at all times may be ob- 
tained 


(12) 


E 
= 

3 

EvarA + 26.) + EwawAw(n6; + 

— (13) 
ErAr + EwAw 

where » = q,/l and fort < t), 6; = 0, and fort > 4, » = 1.0, and 


From Eq. (13) general expressions for the maximum thermal 


-o- 
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em 


Fic. 1. Typical I-section 


stresses in the flange and web may be obtained. 


or = [Er/3(ErAr + EwAw)] |ErAr(@: — &)ar + 
EwAw + — (14) 
+ EwAwaw(n6; — (15) 


Since durirg heating or is compressive and ow is tensile, it 
may be required to investigate thermal buckling of the flange. 


For this, the average flange thermal stress is required. Hence 
from Eq. (13) 
or ave = |ErEwAw/3(ErAr + EwAw)| X 

law(n0; + 263) — + 262)} (16) 


Formulas have been developed for an I-section giving the 
maximum and average thermal stresses in the flange and the 
maximum web thermal stress in terms of structural and material 
properties, and three temperature parameters. These parameters 
can be approximately defined by relatively simple functions also 


presented. 
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Duration of a Constant-Mass Aircraft With a 
Given Airframe Weight and Arbitrary 
Energy-Source Weight 


Henry R. Jex 
Preliminary Design Department, Radioplane Division, 
Northrop Aircraft, Inc., Hawthorne, Calif 


March 28, 1958 


A’ EQUATION for the duration of a constant-mass aircraft 
whose energy is supplied by a fixed-weight source was de- 
rived years ago by technically inclined aeromodellers* for deal- 
ing with the endurance of a model aircraft powered by elastic 
rubber bands. This note reviews that equation and considers 
its implications for nuclear-powered aircraft. 

Assume an aircraft of constant weight (I’), flying at some 
constant speed (J), with a constant propulsive efficiency (). 
It does not matter whether or not there is a climb-plus-descent 
as long as the starting and landing altitudes are equal, the flight 
speed is constant, the path angles are small, and the air density 


* Including W. H. Phillips, E. E. Larrabee, and Alan C. Brown 
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(p) does not vary appreciably. The power is assumed to be 
supplied with average efficiency (7) from a source having arbi- 
trary but constant weight (W,), and having associated with it a 
specific energy-storage capacity (E = total energy available + 
Wy). This “energy-source weight’’ is that portion of the gross 
weight which is proportional to the energy stored—e.g., the 
weight of rubber bands in a model aircraft, or prorated weight 
of the reactor plus auxiliary components in a nuclear-powered 
The remainder of the gross weight is termed the “‘air- 
The total weight (W = We + W’4) is 

We examine the variation of duration 


aircraft. 
frame weight” (W"4). 
constant during flight. 
(7) fora given W4 as is varied. 

The duration is given by the basic relation: duration = 


propulsive energy available/velocity X drag. From this, 


T = E We n/V (W) (1) 
With the additional relationships for flight speed and weight, 
V = [(2/p) (W/S) (1/Cx)|"?, W= Wet+ Wa 
we get (by manipulation) the equation for duration: 
T = nE(p/2)"? [We/(We + Wa)*’?] (2) 


We find that the endurance is proportional to nF as expected, 
and to (C,* ?/Cp), which is the familiar minimum-power param- 
eter. It is interesting to note that for maximum endurance the 
flight should be performed at low altitude and air speed (high 
p; low V) and thus the original assumption of nearly constant 
air density is justified. 

Let us now introduce the ‘“energy-to-airframe-weight ratio” 
(K = We/W4,), and the “bare-airframe wing loading’ (W'4/S), 
and rearrange Eq. (2) as follows: 


T = f(K) (3) 


where 


Eq. (3) gives the effect of adding an increasing percentage of 
energy source to an already existing basic airframe of weight 
Ws. By inspection or differentation, it is found that f(A), and 
hence 7, have a maximum for A = W¢/Wa4 = 2.0 (see Fig. 1). 
That is, when adding a constant-weight energy source to an 
existing airframe, the maximum endurance is obtained when the 
energy source is twice the weight of the bare airframe. The physical 
reason is obvious: for A = O there is no energy available; for 
K = o@ the total wing loading, required velocity, and, hence, 
power required also become infinite so that endurance is_ nil. 
At some intermediate value there must be an optimum. 

Besides being of use to aeromodellers, this result is also of 
interest to the designer of long-endurance nuclear aircraft (radar- 
picket airplanes, for example) where take-offs and landings are a 
minute fraction of the overall endurance and a long time be- 
tween nuclear recharging is desired. The conclusions drawn 
from this equation are: 

1. Flights should be performed at the lowest altitude possible. 

2. Fly at or near the maximum ratio of C,* ?/Cp, which im- 
plies: (a) Flight-path control will be difficult due to the altitude- 
control reversal effect below minimum-power speed. (b) High- 
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employed. 

3. Low-wing loadings are desirable. 

4. The energy-supply weight should be at least 50 per cen; 
but not more than 67 per cent of the gross weight. 

5. Due to cube-square-law effects on the airframe wing 
loading, and exclusive of the power-plant-shielding problem, the 
maximum endurance could probably be attained with moder. 
ately small aircraft. 


Detached Shock Waves Ahead of Gas-Sampling 
Probes 


Robert Friedman* and Donald R. Boldman** 
Lewis Flight Propulsion Laboratory, NACA, Cleveland, Ohio 
March 31, 1958 


SYMBOLS 


C = slgtan — VB? tan? yg, — 1) 

L = distance from detached shock to probe along centerline, in 

Mach Number 

yj = probe inside radius, in 

ry = probe outside radius, in 

Vm = free-stream location of streamline that separates the mass flow 
entering the probe from that passing outside the probe, in 

ys = distance from probe centerline to sonic point on detached shock, jn 

B= VM?-1 

n = angle between sonic line and normal to free-stream flow direction 

¢s = slope of detached shock at sonic point 


F’ DETERMINATION of the approach to chemical equilibrium 
in exhaust nozzles, gas composition is measured in_high- 
energy streams, that is, regions of high temperature and high 
Mach Number. A representative sample is extracted from the 
stream through a cooled probe and then passed to the analyzer 
It is assumed that the gas composition at the analyzer corre- 
sponds to that of the free stream ahead of the probe. If a de- 
tached shock wave forms ahead of the probe, however, the gas 
composition sampled by the probe will no longer be that of the 
free stream. Reaction will continue to take place at greatly in- 
creased static temperatures and pressures in the length between 
the detached shock wave and the gas-sampling probe. The 
reduction of a detached shock to an attached shock would limit 
this reaction length to the shock thickness, or the order of a few 
molecular paths. 

Although it may be possible to design cooled gas-sampling 
probes that will maintain an attached shock wave at all super- 
sonic conditions, a graphical-extrapolation method holds greater 
interest. If shock-detachment distance were proportional to 
probe radius, then extrapolation of data from probes of different 
radii to those corresponding to a zero-radius probe would give 
data effectively taken at zero detachment distance, or an attached- 
shock condition. 

As a first step, the feasibility of the extrapolation of detach- 
ment distances was investigated. Measurements of the shock- 
wave-detachment distances in front of bluff-nosed probes were 
taken. For this purpose, a multiple-probe holder was installed 
in a 610°R. total temperature, Mach 1.9 wind tunnel, and the 
individual probes were connected to a vacuum pump to simulate 
While the tunnel could not simulate 4 
high-energy environment, it observation of gas 
sampling probes optically by an interferometer. The shock- 
detachment distances ahead of the probes were measured from 


gas-sampling operation. 
permitted 


the interferograms such as shown in Fig. 1. 
For correlation of the data, a shock-detachment analysis 4s 
given by Moeckel! was used: 


* Aeronautical Research Scientist. 
** Co-op Student, University of Cincinnati 
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Fic. 1. Interferogram of simulated gas-sampling probe at 
Mach 1.9. 
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Fic. 2. Sketch of detached shock ahead of gas-sampling probe. 
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Fic. 3. Experimental location of detached shock waves. 
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L/r, = (C + tan n) — tan — Cly»,/r, (1) 


The symbols of Eq. (1) are illustrated by the sketch in Fig. 2 


where the shock shape is drawn from the experimental inter- 
ferogram (Fig. 1). The sonic point on the shock y, was located 
from the measured shock slope and conical flow charts?; from this 
point a straight line was drawn to the sonic point on the body, the 
probe shoulder. Although considered a sonic line in Moeckel’s 
analysis, this line was a control line’ which established the values 
of y and did not necessarily coincide with the actual sonic line cal- 
culated from the density ratios determined from the interfero- 
grams. The best value of y,, for substitution into Eq. (1) was 
found to be the inside radius of the internal connecting tube 
within the probe tip (Fig. 2). These small connecting tubes be- 
tween the probes and the vacuum pump were required to ac- 
commodate several probes in a single multiple-probe holder 
Normally the probe inside radius would be used for y,,. 

The experimental points are plotted in Fig. 3. As each probe 
had a different value of y,,/r,, a family of lines with the same inter- 
cept at the origin but varying slope was constructed from Eq 
(1) to fit the data. Each experimental point shows excellent 
agreement with the corresponding correlating line in Fig. 3; in 
fact, the points could be correlated with little scatter by a single 
line. Nevertheless, it would be important for an accurate extra- 
polation in multiple-probe gas-sampling installations to keep the 
probes free of internal constrictions and to have the same r,/r, 
ratio for each probe 

The results thus indicated that, for a given set of bluff-nosed 
probes, extrapolation to an effective attached-shock condition is 
possible. The extension of these results to actual gas composi- 
tion measurements appears encouraging; for in the in- 
compressible, constant-temperature region behind the shock, 
reaction rate is proportional to the detachment distance, and the 
correlation of measured gas composition would resemble that of 
detachment distance. 
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On the Performance of a Double-Diaphragm 
Shock Tube Using the Reflected-Shock 
Method and a Light-Gas Buffer 


Charles J. Schexnayder, Jr. 
Aeronautical Research Engineer, Langley Aeronautical 
Laboratory, NACA, Langley Field, Va. 


April 1, 1958 


iy THE PAST FEW YEARS a need has developed for hypersonic 
shock tubes. Many different methods had previously been 
proposed for generating strong shocks in shock tubes.' One 
combination of these was utilized by the Langley Laboratory of 
NACA in 1956 when it built and put into operation a double- 
diapbragm shock tube to study the effects of high temperature 
on air properties. 

Recently Hertzberg? has proposed interesting configurations 
to extend the usefulness of shock tubes in high-temperature re- 
search. It is hoped that this note may provide some useful in- 
formation regarding the performance of a shock tube of this 
more recent type. 

The shock tube is shown in Fig. 1 and consists of a conventional 
4-in. diameter shock tube (sections I and II) driving a smailer 
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Ist Diaphragm 2nd Diaphragm 


Driven Gas, Air O5 To 5mm) 


Section | m Section Il 


(Py) (Pp) 4°D : 
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Fic. 1. Double-diaphragm shock-tube configuration. 
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Fic. 2. Comparison of theoretical and experimental shock- 
wave strengths for double-diaphragm shock tube. Optimum 
theoretical values of ?/7?; indicated. Helium is used for buffer 
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Fic. 3. Comparison of experimental attenuation of shock waves 


using various shock-tube drivers. 
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l-in. diameter shock tube (section III) by means of the “fe. 
flected-shock”’ method, which essentially allows for the maximum 
shock-tube-driver temperature due to the total reflection of the 
shock wave in section II and the most efficient subsonic expansioy 
of the driver gas in the contraction of section II. Rupture of the 
second diaphragm is caused by the reflection of the first shock 
wave from the end of the buffer section. At this time the pres. 
sure difference across the diaphragm suddenly increases from 
about one atmosphere to several hundred psi. It has been found 
from simultaneous pressure records in the buffer chamber ang 
ionization-gap detector records in the shock tube for the case of g 
helium driver that this rupture time is negligible for (.008-jn 
spring brass scribed into eight pie-shaped sections, (0.004 to 
0.005 in. deep. For the high-pressure diaphragm, steel or stain. 
less steel is used, scribed into four pie-shaped sections. The 
scribing of the metals prevents pieces from being torn off. 

The theoretical performance of this configuration is readily 
determined for the nonviscous case if the contraction ratio at the 
second diaphragm is assumed infinite so that steady-flow expan- 
sion to Mach Number 1.0 may be assumed at the throat as well 
as total reflection of the buffer shock wave. It is to be noted 
that for each given gas combination and overall diaphragm pres. 
sure ratio there is an implicit buffer-chamber pressure for maxyi- 
mum driven shock Mach Number providing no other boundary 
restrictions are imposed on the diaphragms. (Hertzberg, in 
reference 2, fixes the pressure difference across the second dia- 
phragm for mechanical reasons and finds an optimum gas com- 
bination. ) 

The theoretical performance of the double-diaphragm shock 
tube for two gas combinations has been computed and is shown 
in Fig. 2 as the variation of shock Mach Number in section II] 
and as a function of the overall diaphragm pressure ratio. The 
hydrogen-driven case is computed for the case of real hydrogen 
according to reference 3. The experimentally determined maxi- 
mum shock Mach Numbers found from the shock velocity meas- 
urements (see Fig. 1) for the two gas combinations and two dif- 
ferent second diaphragm materials are also shown. 

Shown in Fig. 3 is the variation of shock Mach Number along 
the shock tube (attenuation) for some typical tests, and for com- 
parison, data from references 2, 4, and 5. By comparing the 
attenuation curves for the hydrogen and helium drivers, one can 
see that a considerable increase in shock velocity has been 
achieved for testing purposes by using hydrogen instead of 
helium as a high-pressure driver. 

It should be noted, however, when examining Figs. 2 and 3 that 
the maximum experimental velocities obtained in section III 
for the two gas combinations were lower than theoretically com- 
puted. For the helium driver, no actual maximum velocities 
were obtained so that the highest shock velocities recorded were 
used in Fig. 2. For the hydrogen driver, peak velocities were 
found at a considerable distance (90 diameters) from the second 
diaphragm. This could have been due to a lower buffer gas 
pressure and a lower reflected shock pressure in section I], as 
compared to the helium case, thereby allowing for a slower dia- 
phragm opening time. No pressure records were taken for the 
hydrogen case. 
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Optimization of Multiweb Beams Under 
Combined Bending and Torsional Loading 


Gabor Strasser 


Structural Research Engineer, Special Weapons Division, Bell 


Aircraft Corporation, Buffalo, N.Y. 
April 11, 1958 


SYMBOLS 
Vv applied bending moment 
7 applied twisting moment 
chordwise length 
beam depth 


PTIMUM CONFIGURATIONS and optimum stress-envelopes for 
O multiweb beam construction under pure bending were in- 
vestigated by G. Gerard! and A. Schnitt,? et al., respectively. 
The purpose of this note is to report briefly the findings of an in- 
vestigation dealing with the extension of this work (see Fig. 1). 

Since the ultimate load condition on a given multiweb when 
used in a wing structure may be a combination of bending and 
torsion, it is of interest to determine the effects of combined load- 
ing on Optimum proportions as well as on optimum stress. 

The method of approach used was similar to the ones in ref- 
erences (1) and (2). A twisting moment was added to the bend- 
ing moment however, and this combination constituted the overall 
loading condition. A bending moment was always assumed 
while the twisting moment was considered a variable, going from 
zero to a Value equal to the bending moment. By an iterative 
process, it was shown that the outer two webs transfer at least 83 
per cent of the maximum shear flow between cover skins. The 
remaining 17 per cent is distributed among all the remaining in- 
termediate webs. Therefore, shear loading of the intermediate 
webs may be neglected for all practical purposes in the optimiza- 
tion process, and the following loading conditions can be assumed 
for the plate elements of multiweb beams: 

(a) Skins are loaded in combined compression or tension and 
shear. 

(b) Inner webs are loaded in bending only. 

(c) Outer webs are loaded in combined bending and shear. 

Under these conditions, the requirements placed upon the 
beam were that all elements buckle simultaneously, thus prevent- 
ing the retention of reserve strength after failure*; and the solid- 
ity of the cross section must be a minimum. 

It was found that while the absolute dimensions of the elements 
increase with the addition of the twisting moment, the relative 
optinum proportions change but little from the findings based 
on bending only. In fact, they change so little that for all prac- 
tical purposes the same optimum proportions can be applied to 
the combined loading condition, except in the case of the outer 
webs. (See Table 1.) 

To obtain optimum stresses under combined loading, the envel- 
opes as presented in reference 2 may be used with the redefinition 
of the structural index on the abscissa. But here again the varia- 
tion will be small, as it will never exceed 5 per cent. (See Table 2.) 

It can be concluded, therefore, that the optimum proportions 
of multiweb beams (except the thickness of the two outer webs) 

*The foregoing analysis did not include post-buckling strength which 
was shown to exist by Rosen‘ in the region where the critical stress is less 


than 75 per cent of the yield stress 


UPPER SKIN 


= 
| 
4 
J 
WER SKIN INNER WEBS 


Fic. 1. Cross section of multiweb beam. 


FORUM 529 


TABLE 1 
Tabulated Optimum Proportio: s 


Loading Bending Combined Bending 
Proportions and Torsion 


Shape of cells Approximately Approximately square 


square 
Inner web 40% of skin 38-40% of skin thickness 
thickness thickness 
Outer web 40% of skin 40% of skin thickness when 
thickness thickness T/M = 0, and 90% of skin 


thickness when 7/.\/ = 1 
Straight line interpolation 
may be used for inter 
mediate value of 7/1. 


TABLE 2 
Comparison of Structural Indexes 


Bending Only Bending with Torsion 
M/ch? (M/ch?) {1 — 0.05 


and applicable optimum stress envelopes? are insensitive to tor- 
sional loading. This is true as long as the twisting moment does 
not exceed the value of the bending moment. Since this restric 
tion still encompasses the practical range of loading combina- 
tions as applicable to wing structures, the conclusion may be 
generalized as follows: 

Optimum proportions are independent of material properties 
and the results based on a practical combination of bending and 
torsional loading do not vary appreciably from the results ob- 
tained in reference 1 and 2 where a pure bending moment was 


considered only. 
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Minimum-Drag Cone Frustum at Hypersonic 
Speeds 


Robert W. Truitt 

Professor of Aeronautical Engineering, Virginia Polytechnic 
Institute, Blacksburg, Va. 

April 14, 1958 


—— PAPER GIVES an aerodynamic interpretation of a minimiza- 
tion problem which was first solved geometrically by R. L 
Ellis.!. The results are shown to be of importance in the design 
of hypervelocity missiles. 

Consider a general pointed body of revolution with the fixed 
nose coordinate A(O0, 0) and the fixed shoulder coordinate B(x», 
ye). It is known that, under the assumption of simple impact 
theory, the zero-lift pressure drag on this body is practically the 
same value as that for its inscribed cone. That is, neglecting 
centrifugal forces, any pointed body of revolution in hypersonic 
flow can, for practical purposes, be replaced by its equivalent in- 
scribed cone of semivertex angle 6 (see Fig. 1). Furthermore, it 
can be shown by the calculus of variations that the drag on the 
body AB of given fineness ratio can be reduced by replacing it 
with the cone-frustum body A 7B (see Fig. 1). 

The pressure-drag coefficient on the body AB can be written, 


= 
+ 
= 
High 
pp. | 
I 
ib 
hock § 
TN 2 
sley | 
= 
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y 
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EQUIV. 
tan 
Je tand 


M >> 1 8 

tan 8 | 


Fic. 1. 
BODY RADIAL 
DISTANCE FR=5 
1.0 
8 
Newton 
body 
6 min.-drag 
cone frustum 
4 
%4-power body 
0 


ie) 2 4 6 8 1.0 
BODY AXIAL DISTANCE 


Fic. 2. 


according to simple impact theory, as! 
x2 
Dan = 27 0 Ly dx /( 1+ | 


where y is the local body ordinate, p = dy/dx is the local slope 
of the body, and / is the integral to be used in the minimization 
procedure. 

Now the change in drag in passing from the body AB (or its 
equivalent cone) to the cone-frustum body A 7B not only involves 
changes from J to J + 6/ but also results in the appearance of the 
drag on the normal surface A 7’ = 6y. Thus, the total change in 
drag in going from the body AB to the cone-frustum body A7B 


= 27/1 (1) 


may be represented as 
ACp = 2 yo by) + 2 (2) 
By the usual methods of the calculus of variations, the first-order 
terms in 6/, according to Eq. (1), may be expressed as! 
= — [y(3p? + p*) dy/(1 + — 
[2p? sy dx/(1 + p)2] (3) 


Substituting this in Eq. (2) and simplifying, 
ACp Yo by = po?) + po?)?] 
x2 
f, [2 éy dx/(1 + po?)?] (4) 


where po = tan @ is the slope of the frustum surface 7B. It is 
evident from Eq. (4) that if, for example, pp = 1 (@ = 45°) the 
drag is reduced in passing from the body AB to the cone-frustum 
body ATB. On the other hand, the maximum drag reduction in- 
volves a functional relation between the slope po and the given 


fineness ratio, x2/2y2. The following simple method can be used 


1958 
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to find this functional relation 
6 = 646) 5 
and thus determine the minimum-drag cone-frustum body, 
The total drag coefficient for the cone frustum, Cy,, is (see ) 
Fig. 1) 


Cor = Coar + Corp 


where Cy47 and Cpr, are the drag contributions of the normal ang 
slant surfaces of the frustum. Therefore, based on unit area oj 

the base (7 y.? = 1), 


Cpr = 2[1 — (tan 6/tan 8) |? + 2 sin? 9}1 — [1 — } 


(tan 0/tan 6) (7 


Setting to zero the derivative of Eq. (7) with respect to @ (keeping 
6 = constant) gives the minimum-drag cone frustum described by 
the relation [see Eq. (5) | 


tan 26 = 2tané (8 


The problem of aerodynamic heating dictates the use of a blunt 
nose section of fineness ratio of the order of unity. For example, 
a hemispherical nose has a fineness ratio of 0.5 (6 = 45°); thus, 
by Eq. (8), @ = 31.7° and from Eq. (7) Cog = 0.77. Since the 
drag coefficient for the hemisphere, neglecting centrifugal effects, 
is 1.0, this represents a 23 per cent reduction in drag. Fig, 2 
shows that, for this typical small fineness ratio, the cone frustum 


is a good approximation to the exact minimum-drag shape (New- | 


ton’s body!). Fig. 2 also shows that, for small fineness ratios, 
the 3/4-power body is no longer a valid approximation to the 
minimum-drag shape. 

Eq. (7) gives a simple solution for the drag coefficient which, for 
practical purposes, is that for the exact minimum-drag shape of 
small fineness ratio. Finally, from the point of view of construc- 
tion and heat-transfer reduction, the geometrical simplicity of 
the minimum-drag cone frustum has important design implica- 
tions. 
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Supersonic Flow Near the Junction of Two 
Wedges 


Franklin D. Hains 

Aerodynamicist, Reseerch Division, Bell Aircraft Corporation, 
Buffalo, N.Y. 

April 14, 1958 


pee LINEARIZED SOLUTION for supersonic flow near the inter- 
section of two quarter-infinite wedges was obtained by 
Snow.! The purpose of this note is to describe the exact solu- 
tion to this problem. 

The configuration is shown in Fig. 1 with the x axis parallel 
to the free-stream Mach Number J/,. Also shown is the Y}Z 
coordinate system whose Y axis lies along the intersection of the 
top surfaces of the wedges and whose Y axis lies in the plane of the 
surface of the wedge I. The wedges are assumed to have a plane 
of symmetry which passes through points O, D, R, and C. For 
simplicity, the sweepback angle is taken to be zero, but the 
method of solution can readily be extended to swept leading 
edges. 

A typical streamline is shown in Fig. 2. 
sects shock II at a and emerges parallel to the surface of the wedge 
1. The straight-line projection f’a’b’ is inclined slightly toward 
the Y axis because J, is not parallel to the X axis. At 5 the 
streamline intersects shock III whose orientation is not the same 
as II’ because My <M, but is nevertheless plane because of the 
The flow in-section bc -is uni- 


The streamline inter- 


uniform flow region ahead of it. 
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form and is parallel to the axis of symmetry but not to the wedge 
surface. In order to bring the flow tangent to the wedge surface 
I, the streamline segment cd passes through the Prandtl-Meyer 
expansion ONJE. The flow along de is now tangent to the 
wedge surface I but is directed toward the Y axis. At e the 
streamline intersects the characteristic conoid O/JQ and enters 
, conical-flow region where it gradually turns until it lies in the 
plane of symmetry The streamline is made up of straight-line 
segments except along cd and eg where it is curved. 

In the following analysis it is assumed that 6; < 6 < 7 where 
g, is the value for which points J and J’ coincide and the conical 
pattern breaks down. From /; and 6, the angle of shock II 
and VW» are determined. The component of 1/2 normal to shock 
II is 
M, = KM.jcos 6 tan ¢ + sin 6 [(cos @ + 1) ese 6 — 

tan gcot (1) 


where 
K = }sec? ¢ + [(cos @ + 1) ese 6 — tan cot (2) 
In order for JJ, to be parallel to the plane of symmetry ODC, the 
relation 
KM, }|(cos@ + 1) cse @ — tan ¢ cot X 
(cos + 1) + sin 0} + KM,,} [tan g cot — 1] X 
(cos 6 + 1) esc @ — 1f + Me sin @(cos6 + 1) = (3) 


where .V/;,, is the normal Mach Number behind shock III, must be 
satisfied. The angle ¢ which is determined from these relations 
fixes the orientation of the shocks after they intersect. 

The shocks are followed by a conical expansion region similar 
to that treated by Fowell.2 * isa Prandtl Meyer expansion 
centered about OE. It is convenient to introduce a cylindrical 
coordinate system where L is the axial distance from 0, 7 is the 
radial distance from L and ¢ is the angle r makes with the wedge 
surface /. In the ‘free stream” before the first expansion wave 
OEN where & = a, the velocity components are 


= (a3.M; cos B)/a, = 0, = sin B)/a, (4) 
where ¢ is the velocity potential, a the local speed of sound, and 8 
the angle J, makes with OF. The velocity potential in the fan 
is 
= rV (asMz sin B/a;) + [2/(y — 1)] X 

sin — 1)/(y¥ + 1) la + (7/2) —E+ «| ) + Lo, (5) 


Fic. 1. Flow structure near intersection of wedges. 
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Fic. 2. Typical stream line 
where 
= tan 14/(M; sin B)? — 1+ V (4 + 1)/(y —1) X 
tan + 1)/(y — 1)] [CV sin 8)? — 1] (6) 


After the last expansion wave OEJ, (& = n), the velocity is 
Expressed in the XY, Z co- 


parallel to the wedge surface / 
ordinate system, the velocity potential after the last wave is 


og, = [X cos @+ sin + sin — Vecos (7) 
M, = (9,2 2 (8) 


The velocity components parallel to Y, VY, Z are u, v, w, re- 
spectively. Conical coordinates are introduced, defined as 


VY = Y/X, Z2=2/X (9) 


where 


The equation of the parabolic line RNJQ which divides the hyper- 
bolic and elliptic regions is 
(1 + Y2 + Z2) + v2? + w? — a?) = (uw + + wZ)? (10) 
This line is composed of the three segments RN, NJ, and JQ be- 
cause of differences in the velocity impinging on the line. 

The potential equation is elliptic in the region bounded by the 
parabolic line and the wedge surfaces. A new function f(¥, Z) is 
introduced, defined by the relation 


= XM[1 +f) 1+ mY + nZ] (11) 
where /, m, n, are the direction cosines of the vector ./; with 
respect to the Y, Y, Z axes. The potential equation becomes 
+ (+ — 2VZfrz} = 

+ Gfyy + 2FGfyz (12) 
where 


(fy + m)? — (fz + m)*]/2 (13) 


[Zl (01+ Vfy Zfz—- Ufz+n)) (14) 
The velocity components are known along RNJQ. Along 
OC, w = 0, and along CR the velocity component normal to CR is 
zero. These boundary conditions enable Eq. (12) to be solved 
by the relaxation method. 
REFERENCES 
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Effect of Air Pressure on Vortex-Shedding 
Frequency of Cylinders 


Reynold F. Rimoldi,* William D. Clifford, * and Robert J. Bacigalupi* 
Lewis Flight Propulsion Laboratory, NACA, Cleveland, Ohio 
April 18, 1958 


— VORTEX-SHEDDING FREQUENCY from circular eyvlinders is 
generally considered to be a function of the Reynolds Number 
(Re), the cylinder diameter (d), and free-stream velocity (2). 
Previous investigations have correlated the Strouhal Number 
S = fd/u as a function of the Reynolds Number of eylinders. 
In the past, several investigators (references 1, 2, and 3 for 
example) studied the vortex-shedding phenomenon under stand- 
ard conditions of atmospheric pressure and room temperature. 
Data from these investigations were grouped into three distinct 
Reynolds-Number ranges and an empirical Strouhal-Reynolds 
Number correlation was found by Roshko. In most previous 
research, the free-stream velocity and cylinder diameter have 
been common variables and consistently agreed with Roshko’s 
empirical equations. In a more recent investigation, the effect 
of temperature was found to be not completely accounted for by 
the Reynolds Number.! The major purpose of this investigation 
was to determine whether Reynolds Number completely ac- 
counted for the effect of air pressure on Strouhal Number. 

Air from the central laboratory system was passed through an 
electric heater and metered through a standard ASME orifice. 
The air was then passed through a calming chamber and into a 
3- by 45/s-in. test section. Pressures from 4 to 40 Ibs. psia were 
controlled by an exhaust system at the downstream end of this 
test section. A 0.0625-in. diameter rod was placed in a trans- 
verse horizontal position at the center of the test section, and 
shedding frequencies from this rod were sensed by a constant- 
temperature hot-wire anemometer located at a position 5 diam- 
The electrical signal passed through a high- 
pass low-pass filter and an amplifier then fed to the vertical axis 
Frequency was deter- 


eters downstream. 


of an oscilloscope and to a wave analyzer. 
mined by connecting an oscillator to the horizontal axis of the 
oscilloscope so that the vortex-shedding frequency could be deter- 
mined by the Lissajous-figure method and checked with the wave 
analyzer. 

Within the ranges of this investigation, Reynolds Number 175 
to 1,000, and pressure 10 in. to 50 in. of mercury absolute, a defi- 
nite variation in Strouhal Number was observed. At constant 
Reynolds Number, the Strouhal Number rose with increasing 
pressure from lower values and approached the values observed by 
previous investigators. However, as the pressure was further 
increased to 50 in. of mercury absolute, the Strouhal Number 
began to decrease in value. Results of this investigation show a 
definite dependence of the Strouhal Number of cylinders on pres- 
sure and suggest the need for a new correlation including the 
effects of temperature and pressure in addition to the effect of 
Reynolds Number. 


* Co-op Student, University of Detroit. 
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The Influence of a Local Cross-Section ty - 4 
Change on the Two-Dimensional 
Wave by Linearized Theory La 


Ingolf Teipel 
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with th 
. FROM the linearized differential equation for one. a 
dimensional unsteady flow, the influence of locally increasing | aa 

or decreasing cross section on a plane pressure wave will be calcu- 
lated. Three different domains will be distinguished (Fig. 1); | where 


(1) The domain of constant cross section, x < x. the loc: 
(2) The domain of conical cross section, x; <x < x. x 


(3) The domain of constant cross section, x2. < x. At x 


The following relations hold for the velocity potential ¢ in the inasin 
three different domains: The 
sented 

= F(x — at) + F(x + at) 

= (1/x) — act) + (1/x) + act)? (1) ginning 

= F(x — act) + Fe(x + aot) 

a certa 


where the F, are arbitrary functions of their respective arguments | 
compa 


and 
to be e 
= U, —(l/ander = [2/(y — 1)] (a — aw) Asc 
The boundary conditions for F; and F, ia terms of F, and F, are been ¢ 
atx =: having 
by abc 


U(x1,t) = Fi'(x1 — at) + + at) = sntervs 
—(1/x1?) — aot) + at) — 
(1/21?) + aot) + + aot)) (9 

[2/(y — 1)] — ao] = — aot) — + 


Becaus 


= — aot) — (1/1) + aot) mainin 
where the primes denote the differentiation with respect to the eas 
arguments of the functions. The disturbance of the local sound seein 
speed represents in the linearized theory the pressure disturb- 
ance. At x = x. corresponding relations are inserted. 
A saw-tooth simple pressure wave (Fig. 2) is moving into the ‘ioe 
fluid at rest, supposed to be the initial condition (f = 4): New Y¢ 
x — at) = | 
A[(x — xo)/b], + act) =07 (3) | 
x< xX, xe: B(x — at) =0, + at) =0 
A is the amplitude of the wave. At x = x, Eq. (2) has to be 
satisfied by the general form of Eq. (3), which holds for the time 
to < t ty + 
F,'(x — aot) = Al(x — aot — xo + (4 


Because the fluid ahead of the wave is at rest, the function F, is 
equal to zero. Buta reflected plane wave F,arisesatx = | 
to <t Sto + (h/ao): — aot) + 
Fo'(x1 + aot) = —(1/x1?)F3(x1 — aot) + 
(1/x1)F3"(x1 — aot) 
F,'(x1 — aot) — + act) = (1/x1)F3’(x1 — aot) 


‘ Fic. 3 
The solutions of Eq. (5) with the initial conditions Eq. (3) are: 
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Fic. 1. Increasing cross section. 


— SX — aot LM — aslo: 


F3(x — aot) = [1 + exp [(1/2x1) (x — aot — x1 + aolo)] — [(x — aot — x0 + — (6) 


+ aol < x + dol li + adolo: 


+ aot) = —A}[1 + 2(x1/lo)] exp [—(1/2m1) (x + aot — x1 — aoto)| — 40x1/lo) + 2(x0/lo) + [C(x + aot — x9 — 


Fort > ty) + (1y)/ay), Fi’ = O, and Eq. (5) takes on the following form: 


Fo'(x1 + aot) = —(1/x1?) — aot) + — aot), — + aot) = — aot) (7 
with the solutions 
ty — doly > — aot: Fy(x — aot) = [1 + exp [—(lo/2x1)| — 2(x1/lo)} exp [(1/2x1) (x — aot — x1 + aol 
+ lo + auto <x + aot: + act) = —Af[1 4+ 2(x1/lo)] exp [—(lo/2x1)] — 2(x1/lo)} exp [—(1/2x1) (x + aot — xo + aoto) 


where use has been made of the conditions that the velocity U’ and 
the local sound speed are continuous at ¢ = f) + (/)/ay) and x = 

At x = x» a second reflection takes place which can be treated 
ina similar manner. 

The pressure (a — a) and velocity distributions are pre- 
sented in Fig. 3 for / = tf) and t = t) + (4/)/ao).. The amplitude 
A is substituted by the initial velocity U; atx = x. At the be- 
ginning both distributions have the same structure; at ¢ = fy + 
4],/ay) a long drawn depression wave has been developed. In 
a certain domain the velocity assumes negative values. For a 
comparison of the waves in Fig. 3 the respective cross section has 
to be considered. 

As control of the disturbances, the equation of continuity has 
been checked. It is found that the mass within the saw-tooth, 
having constant wavelength ia a linearized theory, has increased 
by about 85 per cent of the initial disturbed mass 1/;. In the 
interval —2.0 < x < + 2.0 the wave contains 60 per ceut of 
JJ, as depression wave. 

1.85.1, — 0.601, — 614 = M, 
Because the perturbation mass must be equal to all times, the re- 
maining part 6,.1/; for the complete compensation is distributed 
in the inner domain. 

From this calculation it can be seen that a long drawn wave 
develops if a saw-tooth wave runs through a variable cross section. 
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i Fic. 3. Pressure and velocity distributions for the increasing 
cross section as shown below. 
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The Effect of Gas Properties on the Heat 
Transfer in Stagnation Flows 


Ivan E. Beckwith 
Langley Aeronautical Laboratory, NACA, Langley Field, Va 
April 21, 1958 


SYMBOLS 


1 + — 


A = , Yaw parameter 
1 + — 1)/2]M.,2 cos? A 
8 = 2m/(m + 1), pressure gradient parameter 
m = (X/U,) (dU-/dX) 
a = transformed’ chordwise coordinate 
l = transformed’ chordwise velocity 
M. = stream Mach Number 
A = yaw angle 
¥ = ratio of specific heats 
S = Sutherland's constant 
7 = temperature 
Subscripts 
9 = total-stagnation conditions 
e = local values external to boundary layer 


r A RECENT sTuDY'! of real-gas effects on the heat transfer at 
the stagnation point of a body of revolution, it was concluded 
that the major deviation in the heat-transfer parameter from the 
low-temperature, perfect-gas value is due to the variation of the 
product of density and viscosity py across the boundary layer 
In particular, for equilibrium dissociation and Lewis Number 1, 
for which the diffusion terms disappear from the differential 
equations,' the heat-transfer parameter was given by the simple 
correlation equation 


Nu V Re = 0.67( petts/ pct (1) 


where the symbols are those of reference 1 

The purpose of this note is to report that Eq. (1) also correlates 
the results of perfect-gas solutions provided that Sutherland’s 
formula for viscosity is used. Furthermore, the heat-transfer 
parameter at the stagnation line of an infinite cylinder in yaw is 
correlated by a similar equation for a large range of stream Mach 
Numbers and yaw angles. 

The boundary-layer equations solved? are the same as those of 
reference 3 except that pu is allowed to vary across the boundary 
laver according to the perfect-gas law and Sutherland's viscosity 
formula rather than py = constant.’ Other assumptions are 
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Fic. 1. Correlation of heat-transfer parameter at the stagna- 

tion point of a body of revolution and at the stagnation line of a 

yawed infinite cylinder. Open symbols are for \ = 1, vertical 

line through symbols denotes \ = 3.0, and a crossed symbol 
denotes A = 11.0. 


constant specific heats, constant Prandtl Number, isothermal 
wall, and similar boundary-layer profiles. The latter assump- 
tion—that is, that the dimensionless velocity and thermal pro- 
files are functions of a single variable—is physically correct at 
the stagnation line of the cylinder. 

Values of the heat-transfer parameter obtained from these 
solutions are plotted in Fig. 1 versus the ratio (p.u./p.u,.) where 
subscript s now denotes values in the external flow at the stagna- 
tion line of the cylinder as well as the stagnation point of a body 
of revolution. 

The solutions for 8 = 0.5 and A = 1.0 are applied to a body of 
revolution by the use of Mangler’s transformation relating two- 
dimensional and axisymmetric flows. 

Since the Sutherland viscosity formula and constant Prandtl 
Number are used in reference 1, as well as in the present solutions, 
the close agreement between the present results for 8 = 0.5 and 
2q. (1) indicates that the heat-transfer parameter at a stagnation 
point is not sensitive to the effects of dissociation on density and 
specific heat within the boundary layer. For equilibrium disso- 
ciation and Lewis Number 1, the heat-transfer rate at a three- 
dimensional stagnation point can then be calculated from the 
equation! 


du = V de/dx), — Ie)/o] (Nu/V Re) (2) 


where the symbols are those of reference 1, and all quantities are 
evaluated in the real-gas flow except (Nu/+/Re) which may be 
obtained from either Eq. (1) or a solution of the boundary-layer 
equations for a perfect gas with Sutherland’s viscosity formula 
provided that the value of p,u./pwue applying to the solution is 
equal to the real-gas value. 

The heat transfer at the stagnation line of a yawed cylinder is 
also calculated from Eq. (2) except that h, should be replaced by 
the adiabatic wall enthalpy given by the equation? 

h, = Naw = ho — (2-2/2) (1 — 1) (3) 
where 7, is the external spanwise velocity and r ~ +/o%* is the 
recovery factor. The values of the heat-transfer parameter on a 
yawed cylinder from the present perfect gas solutions are corre- 
lated (Fig. 1) by the equation 


Nu/V Re = 0.5( pstts/pwpw)™ (4) 
which is within about 4 per cent of all solutions for 8 = 1 except 
those for X = 11.0 and S/T) = 0.2. These latter values are not 


representative of flight conditions since they would occur at large 
Mach Number and large vaw angle with 7) ~ 1,000°R. By 
analogy with the above results for a three-dimensional stagna- 
tion point, the heat transfer at the stagnation line of a yawed 
eylinder in a real-gas flow (with equilibrium dissociation and 
Lewis Number 1) may be calculated from Eqs. (2), (3), and (4) 
where all quantities are, of course, evaluated for the real-gas 
conditions. 

Perfect-gas solutions for the three-dimensional stagnation point 
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have also been given by Mark,‘ who used the Sutherland viscosit, 
formula. It can be shown that Mark's results are also correlate; 
by Eq. (1). 
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Poisson’s Ratio for Honeycomb Sandwich 
Cores 


George A. Hoffman 
Aeronautics Department, The RAND Corporation, 
Santa Monica, Calif. 


April 28, 1958 


A’ INTERESTING mechanical property of honeycomb cores js 

their Poisson’s ratio: it is very sensitive to cell geometry 
and can assume values from zero to about three. Large values 
may be observed when flexing some slabs of honeycomb with 
flat cell walls, while rippled-wall honeycombs demonstrate zero 
Poisson's ratio. The ratio of anticlastic curvatures is indicative 
of the value of Poisson’s ratio for axial loading in the plane of the 
slab. 

In reference 1 the value of 0.806 is derived for a circular-cell 
honeycomb, while reference 2 obtains a value of almost 1 for 
Poisson’s ratio for regular hexagonal cells. In this note, Poisson's 
ratio is evaluated for honeycomb cores of various cell shapes and 
deformation modes. It will also be shown that large variations in 
the ratio do not affect the stiffness of the sandwich panel.’ 

Equilateral hexagonal and square cells will be considered here; 
the doubled foil line will be the x-direction, and the cell width, ¢, 
will be the spacing between doubled foil lines (see insert, Fig, 1). 
The angle a between the x-axis and an adjacent wall is 60° for 
regular hexagonal cells, though in practice it may be much less 
(the so-called under-expanded core), or well beyond 60° (over- 
expanded cores). Reducing the doubled foil line to a minimum 
results in a nearly quadrilateral cell herein assumed square 
Conforming to common practice, the doubled foil line will be 
assumed parallel or perpendicular to the loading or bend axis. 
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Fic. 1. Poisson’s ratio for hexagonal cell honeycomb core. 
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For uniaxial stress in the middle plane of an unrestrained core 


plate, Poisson’s ratio may be defined as 


magnitude of strain perpendicular to stress direction 


magnitude of strain in the direction of stress 


Therefore, v, will denote the ratio — (e’,/e’,) fora negative e,’ re- 
sulting from y-wise tension stressing, while vy, = —(e,"/e:") for 
negative ¢€,” for x-wise stressing. The derivation of v, and v, will 
be based on uniaxial x- or y-wise stressing, resulting in small elas- 
tic deformations. Bending of the walls of a unit depth core will 
be analyzed, disregarding secondary affects where permissible 
when @ is not near 0° or 90°). A core wall will be considered an 
S-shaped element, a/2 high, composed of the flat portion of the 
foil and of circular arcs at each end, of bend radius r. Since 
the midpoint of the S-shaped element is an inflection point, only 
half of such an element need be considered. 

Applying hypothetical loads P, and P, at this midpoint A, re- 
sults in a moment at any point in the element equal to 


M = P,yx — Pry (1) 
The strain energy, lL’, of the half-element is 


B B 
‘= | (M?2/2EI)ds = (1/2EI (Py%x? — 2P:Pyxy + 


P,2y?) ds (2) 


~ 


Displacements 6, and 6, of the inflection point A, are 


B B 
= 0U/0P, = (Pz En f yds — (P, ED { xy ds 


B B 
= 0U/OP, = (Py ED) x?ds — (P; ED xy ds 


The integrals in Eq. (3) happen to be, respectively, the mo- 
ments of inertia about the x and the y axis and the product of 
inertia about the x and y axes of a line shaped as the half-element. 


> 


(3) 


> 
< 


Neglecting terms in ar? and r* gives 


B 
y? ds = (1/3 sin a) (a/4)3 — (a/4)? [1 — cos + 


B B 
sin a] — af = (1/cot? a) = (1/cot a) xyds (4) 


Then the relation between displacements of A due to P; or Py 
only, simply becomes 


6, = —6, cota (9D) 


The x and y deformations of a periodically repetitive portion of 
the honeycomb core permit finally to derive Poisson’s ratio. It 
can be shown’ that the smallest period of a hexagonal cell is a/2 
(cot a ese aw) long and of width a, while for a square cell it is a 


square, also of width a. With this consideration 


change in period width 


period width 
period length 
, load only) 
change in period length ; 
(6) 
change in period length 
y= 


period length 


period width 
: : ——, (P, load only) 
change in period width 


Using these last definitions and Eq. (5) we obtain 


vy, = cot a (cot a + esc a) = 
1/v, for a hexagonal cell core | (7) 
vy = vy = | for a square cell loaded diagonally 


Expressions (7) are unaffected in the region 70° > a@ > 50° by 
axial wall deformation or by small deviations from straight lines 
and bend ares. In the region 0° < a < 20°, vy reaches a maxi- 
mum and then drops rapidly as shown in Fig. 1 for a particular 
imperfect configuration. As a approaches 90°, terms in ar? and r3 
become significant in the evaluation of »,, to the extent that at 


a = 90° it can be shown? that 
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of sandwich panels. 


vz = OA(a/r) — 0.2 (for fully expanded core (8) 


Fig. 1 shows v, for a > 80° for the particular case of (a/r) = 8 

An analogy may be drawn between a bent core slab and a homo- 
geneous isotropic plate in bending: the ratio of curvatures in 
the x-z and y-s planes is v, for bending with plane stress in the 
x-z plane and is for plane stress. Therefore, Poisson's ratio 
for bending can be taken to be the same as for axial loading, with 
the moment axis, whether x or y, defining which value to use, 
Vz OT Vy. 

It is interesting to show that the large variations in core Pois- 
son’s ratio (hereafter denoted simply as vy.) shown in Fig. 1 have 
negligible effect on a composite sandwich panel. The effective 
Poisson’s ratio for a whole sandwich panel, v,., was derived in 
reference 4 to be 

ve = (yD + vd)/(D+d (9) 


where vy = v of the face material, D and d are the bending stiff- 
nesses per unit width of the spaced faces and of the core, respec- 
tively, in plane strain. 

The stiffness of composite sandwich panels is proportional to 
1/(1 — »,*?) and from Eq. (9) the ratio (panel stiffness with non- 
zero v,/panel stiffness with », = 0) equals 

— [»D/(D + — + vd)/(D + (10) 
giving {[(10) — 1] X 100 as the per cent increase in panel stiff- 
ness when a core with nonzero », is substituted for a core with 
This per cent increase is plotted in Fig. 2 against », for 
various values of (d/D), and for » = 0.33. 

It can be seen from Fig. 2 that large variations of »,, do not 
affect appreciably the stiffness of the panel, since practical de- 


signs seldom go over d/D of 0.001. 


ve = 0. 
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Addendum—‘‘On Axisymmetric Thermal 
Stresses in Thin Shells of Revolution’’ 


G. D. Galletly 
Mechanical and Electrical Engineering Department, 
Shell Development Company, Emeryville, Calif. 


April 23, 1958 


| THE READERS’ Forum, March, 1958 issue, a note appeared 
by the writer in which it was stated that there was an apparent 
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error in an equation given in the book Warmespannungen by 
Melan and Parkus. The writer ascribed this error to the com- 
patibility equation used by them. 

Recent correspondence between the writer and Professor 
Parkus has established that the equations given by Melan and 
Parkus are correct. Those given by the writer are also correct. 
The difference in the two sets of equations is due to the definition 
of the variable Q. Following Timoshenko, the writer used 
Q = rg where qg is the transverse shearing foree. Melan and 
Parkus, on the other hand, used Q = 4roq/6? where 6 is the shell 
thickness. 


Use of the Hydraulic Analogy for ‘‘Inside’’ 
Problems 


R. A. A. Bryant 


Senior Lecturer in Mechanical Engineering, 
The New South Wales University of Technology, Kensington, 
N. S. W., Australia 


May 2, 1958 


HE NOTE! recently published by Syégo Matsunaga provides 

another interesting application of the hydraulic analogy and 
brings to mind the fact that there are a number of steady flow 
“inside” problems which can be conveniently investigated using 
this method. For example, Professor Matsunaga also has used 
the analogy to study the working cycle of safety valves. 

Of course, more pertinent problems such as flow in a variable 
geometry ram-jet diffuser, which necessitates the use of the 
“towed model” analogy, are more difficult to handle,? and it is 
doubtful if quantitative data of known validity will be forth- 
coming from such studies. The problem of transferring unsteady 


hydraulic data for ratio of specific heats, y = 2, to equivalent 
unsteady conditions in air flow where 7 = 1.4, appears to be the 


main difficulty in this type of problem. 

For “inside”? problems it may be noted that two types of 
analogy can be constructed. Of them, a one-dimensional 
analogy, which requires a nonrectangular channel and which 
gives direct predictions for any value of y, is limited to applica- 
tions of an essentially one-dimensional nature but is valid under 
both steady and unsteady conditions. Periodic behavior in engine 
exhaust systems is an example of the type of study which can be 
successfully made with a one-dimensional analogy. 

As soon as the nature of the flow no longer complies with the 
requirements of one-dimensionality, use of a rectangular tank 
cross section is necessary. Consequently studies such as that re- 
ported in reference 1 in which the motion is essentially two- 
dimensional give predictions for the fictitious gas with 7 = 2. 
The validity of such predictions is always in doubt and correla- 
tion with prototype (axisymmetric) behavior is necessary before 


and amplify the aerodynamic loadings by a large 
amount. 

The results also serve to indicate the complex 
manner in which the strains depend upon both the air- 
plane rigid-body and elastic-body motions. The results 
indicate that the reliability of strain calculations will 
depend heavily upon the accurate representation of 
both the magnitude and phase of the rigid-body and 
elastic-body motions. 
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Some Structural Response Characteristics . . . (Continued from page 521) 


extrapolation of hydraulic data can be safely considered. — The 
factors which are most important and require special attentigg 
(other than the inevitable problem of the fictitious gas) may be 
classified under headings of (1) boundary conditions, (2) shock 
conditions, (3) vorticity conditions, and (4) hydraulic flow cop. 


ditions. 

(1) Boundary Conditions——-Boundary-layer conditions age 
never compatible between a physical analogy and its physical 
prototype. This is natural as the motions involved are physically 
different. However, the problem of boundary conditions is more 
insidious as the basic analogy, itself, applies to potential flows 
only. Strictly speaking, we then use a real (water) flow to explaig 
a potential gas flow and there can be no similarity at all with 
respect to boundary conditions. Thus no information aboyt 
boundary-layer shock interaction can be expected. 

(2) Shocks--Shocks break down the assumption of isentropie 
flow which is involved in deriving the basic analogical relation. 
ships. Furthermore, there is no physical correspondence betwee, 
a free-surface discontinuity in water and a compression shock ig 
agas. Normally applications of the analogy to “‘inside’’ problems 
which involve strong shocks vield qualitative data only. 

(3) Vorticity Conditions —It is easy to show that the basic two- 
dimensional analogy extends to rotational motion.* In problems 
where the shocks are straight or slightly curved, it is logical to 
assume a priori that the vorticity distribution will be the same 
(under ideal hydraulic conditions of uniform horizontal velocity 
over the whole water depth) in both analogy and_ prototype 
(potential-flow) systems providing there is no undue influence 
brought about by the side-wall boundary layers in the water flow, 
Vortex sheets (see Fig. 2 of reference 1), however, cannot be con- 
sidered as corresponding since phenomena are involved in the 
water flow which have no place in the analogous gas flow. 

(4) Hydraulic Flow——The hydraulic flow conditions involved in 
the stationary-model technique are a source of partial breakdown 
of the basic analogy mainly due to the boundary layer which forms 
along the bottom of the tank. To overcome this, larger water 
depths are required than that dictated fer free-surface wave 
propagation with zero dispersion. The reader will find further 
information concerning this aspect of the analogy in the paper by 


Itava and Tomita! 
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